In [1]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib

In [2]:
import matplotlib 
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from scipy import stats
from matplotlib.font_manager import FontProperties
from fractions import Fraction
In [3]:
%cd /Users/burbol2/Dropbox/Apps/Computable/Final_Results_CAngles/AnalyticalResults
#%cd /home/eixeres/Dropbox/Apps/Computable/Final_Results_CAngles
/Users/burbol2/Dropbox/Apps/Computable/Final_Results_CAngles/AnalyticalResults

In [4]:
percentages = [0, 5, 11, 17, 33, 50,66]
Waters=[1000, 2000, 3000, 4000, 5000, 6500, 7000, 8000, 9000, 10000]

mac_theta_w = np.zeros(len(percentages))
sigma_w = np.zeros(len(percentages))

mac_theta_m = np.zeros(len(percentages))
sigma_m = np.zeros(len(percentages))

mac_theta_s = np.zeros(len(percentages))
sigma_s = np.zeros(len(percentages))

# In the following dictionaries we will save the raw input data
angles_w = {}
radii_w = {}

angles_m = {}
radii_m = {}

angles_s = {}
radii_s = {}

# In the following dictionaries we will save the results for the cos(theta) and 1/r_base, and their errorbars. 
# We will call them the "results-dictionaries"
theta_w = {}
rbase_w = {}
errortheta_w = {}
error_rbase_w = {}
start_theta_w = {}
end_theta_w = {}
start_rbase_w = {}
end_rbase_w = {}
equil_theta_w = {}
equil_rbase_w = {}

theta_m = {}
rbase_m = {}
errortheta_m = {}
error_rbase_m = {}
start_theta_m = {}
end_theta_m = {}
start_rbase_m = {}
end_rbase_m = {}
equil_theta_m = {}
equil_rbase_m = {}

theta_s = {}
rbase_s = {}
errortheta_s = {}
error_rbase_s = {}
start_theta_s = {}
end_theta_s = {}
start_rbase_s = {}
end_rbase_s = {}
equil_theta_s = {}
equil_rbase_s = {}

# For each surface we will store 10 results corresponding to the 10 droplets in each system. This means that in every 
# entry of the "results-dictionaries" we will store a list of 10 elements. Thus, now we store in them lists with 10 zeros. 

for system in percentages:
    
    theta_w[system] = []
    rbase_w[system] = []
    errortheta_w[system] = []
    error_rbase_w[system] = []
    start_theta_w[system] = []
    end_theta_w[system] = []
    start_rbase_w[system] = []
    end_rbase_w[system] = []
    equil_theta_w[system] = []
    equil_rbase_w[system] = []
    
    theta_m[system] = []
    rbase_m[system] = []
    errortheta_m[system] = []
    error_rbase_m[system] = []
    start_theta_m[system] = []
    end_theta_m[system] = []
    start_rbase_m[system] = []
    end_rbase_m[system] = []
    equil_theta_m[system] = []
    equil_rbase_m[system] = []
    
    theta_s[system] = []
    rbase_s[system] = []
    errortheta_s[system] = []
    error_rbase_s[system] = []
    start_theta_s[system] = []
    end_theta_s[system] = []
    start_rbase_s[system] = []
    end_rbase_s[system] = []
    equil_theta_s[system] = []
    equil_rbase_s[system] = []
In [5]:
# We create an array named "t" with the time values. We use the length of the longest simulation.

maxlength=100 #maximum length in ns
t = np.zeros(2*maxlength)
i=0
for d in frange(0, maxlength, 0.5,closed=0):
    t[i] = d
    i = i + 1
In [6]:
# We save the data in the arrays "theta_all" and "rbase_all"
SAMs=[0,5,11,17,33] # Here we use only the SAMs with Contact Angles. Not the systems with complete wetting.

# Here we create two dictionaries named "angles" and "radii". Each dictionary consists of a set of "key: value" pairs, each one corresponding to a different system. The keys are (also) a pair of numbers "b, c".
# "b" refers to the OH-coverage percentage of the SAM, and "c" to the number of water molecules of the droplet. For example, the key "5, 1000" refers to the system with a surface with 5% OH-coverage, 
# and a droplet with 1000 water molecules. The value corresponding to each key is an array with the 40 calculated angles/base radii, that belong to that system. For example: "radii[(11, 3000)]" would give back
# the array with the 40 calculated values of the base radii corresponding to the system with a surface with 11% OH-coverage, and a droplet with 3000 water molecules.

for b in SAMs:
    theta_data_w, rbase_data_w = np.loadtxt('Contact_Angles2_WaterPeak_s'+str(b)+'.txt', skiprows=2, usecols = (0,1),unpack=True)
    theta_data_m, rbase_data_m = np.loadtxt('Analyt_Contact_Angles2_GDS_s'+str(b)+'.txt', skiprows=2, usecols = (0,1),unpack=True)
    theta_data_s, rbase_data_s = np.loadtxt('Analyt_Contact_Angles2_SAMPeak_s'+str(b)+'.txt', skiprows=2, usecols = (0,1),unpack=True)
    k = 0
    for c in Waters:
        # First data with z=0 at the first water peak
        z = b, c
        th_w = [0] *(2*maxlength)
        r_w = [0] *(2*maxlength)
        th_m = [0] *(2*maxlength)
        r_m = [0] *(2*maxlength)
        th_s = [0] *(2*maxlength)
        r_s = [0] *(2*maxlength)
        i=k*(2*maxlength)
        j=(k+1)*(2*maxlength)
        #print "z=", z, "i=", i, "j=", j
        m = 0
        for l in range(i,j):
            th_w[m]=theta_data_w[l]
            r_w[m]=rbase_data_w[l]
            m = m +1
        angles_w[z] = th_w
        radii_w[z] = r_w
        
        m = 0
        for l in range(i,j):
            th_m[m]=theta_data_m[l]
            r_m[m]=rbase_data_m[l]
            m = m +1
        angles_m[z] = th_m
        radii_m[z] = r_m
        
        m = 0
        for l in range(i,j):
            th_s[m]=theta_data_s[l]
            r_s[m]=rbase_data_s[l]
            m = m +1
        angles_s[z] = th_s
        radii_s[z] = r_s

        k=k+1
In [7]:
z=33,8000
w=11
print angles_w[z][w]
print angles_m[z][w]
print angles_s[z][w]
print radii_w[z][w]
print radii_m[z][w]
print radii_s[z][w]
103.165494112
104.625679364
105.432120474
4.06434237108
4.03879654523
4.02356242953

In [8]:
# func is the function of a line, that will be used in the linear fits
def func(x, a, b):
    return a*x + b

#endpoint returns the integer "end": the last data point different then "nan"
def endpoint(theta):
	for end in range(20, len(theta)):
		 if math.isnan(theta[end]):
			break
	return end

# The following functions will be used to calculate the block averages and the errorbars
# The function naive_variance will only be used inside the function blockAverage
def naive_variance(data):
    n = 0
    Sum = 0
    Sum_sqr = 0
 
    for x in data:
        n = n + 1
        Sum = Sum + x
        Sum_sqr = Sum_sqr + x*x
 
    variance = (Sum_sqr - (Sum*Sum)/n)/(n - 1)
    return variance

def blockAverage(datastream, Nblocks):
    
    # FIRST WE CREATE AN ARRAY TO STORE THE MEAN VALUES OF THE DATA BLOCKS (blockMean)
    blockMean = np.zeros(Nblocks)  
        
    # Nobs is the number of points (observations) in our data
    Nobs = len(datastream) 
    # BlockSize is the size of each block of data
    BlockSize = int(Nobs//Nblocks) 
    #print "BlockSize=", BlockSize
    
    if Nblocks==1:
        errorbar = naive_variance(datastream)
        return errorbar

    else:
       # WE CALCULATE IN A LOOP THE MEAN VALUES OF EACH DATA BLOCK (blockMean)
        for i in range(0,Nblocks-1):
            ibeg = i*BlockSize
            iend = (i+1)*BlockSize
            blockMean[i] = mean(datastream[ibeg:iend])
        # WE TREAT THE LAST BLOCK SEPARATELY, BECAUSE WE HAVE TO TAKE INTO ACCOUNT THE POSSIBLE REMAINING POINTS 
        # WHEN THE NUMBER OF DATA POINTS ISN'T A MULTIPLE OF THE NUMBER OF BLOCKS
        #print "Last block"
        ibeg = (Nblocks-1)*BlockSize
        iend = Nobs
        blockMean[Nblocks-1] = mean(datastream[ibeg:iend])
     
        errorbar = (np.std(blockMean))/sqrt(Nblocks -1) #np.std(blockMean) is the standard deviation of blockMean
        simulavr = mean(blockMean)
        return simulavr, errorbar

#best_start SEARCHS FOR THE STARTING POINT (start) OF BIGGEST TIME INTERVAL where the error interval is smaller then the variation of the data (controlled 
# with a line fit of the data
def best_start(theta,t,omitstart,fixend):

	lastnumber = endpoint(theta)
	endblocks = lastnumber - fixend

	error = np.zeros(endblocks-omitstart) #error will be the total error (of all the points taken each time we choose a different set of data to make the blocks)
	average = np.zeros(endblocks-omitstart) # average is the same, but for the average
	
	slope = np.zeros(endblocks-omitstart)
	intercept = np.zeros(endblocks-omitstart)
	shift = np.zeros(endblocks-omitstart)
	goodblocks = []

	error = error.tolist()
	j=0
	# Loop for finding the "best" interval to do the block averaging
	for i in range (omitstart, endblocks):
		average[j],error[j] = blockAverage((theta[i:lastnumber]),3)
		slope[j], intercept[j], delete1, delete2, delete3 = stats.linregress(t[i:lastnumber],theta[i:lastnumber])
		shift[j] = abs(func(t[i], slope[j], intercept[j]) - func(t[lastnumber], slope[j], intercept[j]))
		if shift[j] <= (2*error[j]):
			goodblocks.append(i)
			#print "i=",i, "shift=",shift[j], "error=",2*error[j]
		j = j+1	 
	if goodblocks==[]:
		#print "not equilibrated yet"
		return None
	else:
		#start = min(goodblocks)+ omitstart
		print "num of good intervals=", len(goodblocks)
		start = min(goodblocks)
		return start
In [9]:
'''
For adjusting the subplots:
subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=None)

left  = 0.125  # the left side of the subplots of the figure
right = 0.9    # the right side of the subplots of the figure
bottom = 0.1   # the bottom of the subplots of the figure
top = 0.9      # the top of the subplots of the figure
wspace = 0.2   # the amount of width reserved for blank space between subplots
hspace = 0.2   # the amount of height reserved for white space between subplots

********************************************************************
For adjusting the legend box:

Padding and spacing between various elements use following
keywords parameters. These values are measure in font-size
units. E.g., a fontsize of 10 points and a handlelength=5
implies a handlelength of 50 points.  Values from rcParams
will be used if None.

=====================================================================
Keyword       |    Description
=====================================================================
borderpad          the fractional whitespace inside the legend border
labelspacing       the vertical space between the legend entries
handlelength       the length of the legend handles
handletextpad      the pad between the legend handle and text
borderaxespad      the pad between the axes and legend border
columnspacing      the spacing between columns

'''
Out[9]:
'\nFor adjusting the subplots:\nsubplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=None)\n\nleft  = 0.125  # the left side of the subplots of the figure\nright = 0.9    # the right side of the subplots of the figure\nbottom = 0.1   # the bottom of the subplots of the figure\ntop = 0.9      # the top of the subplots of the figure\nwspace = 0.2   # the amount of width reserved for blank space between subplots\nhspace = 0.2   # the amount of height reserved for white space between subplots\n\n********************************************************************\nFor adjusting the legend box:\n\nPadding and spacing between various elements use following\nkeywords parameters. These values are measure in font-size\nunits. E.g., a fontsize of 10 points and a handlelength=5\nimplies a handlelength of 50 points.  Values from rcParams\nwill be used if None.\n\n=====================================================================\nKeyword       |    Description\n=====================================================================\nborderpad          the fractional whitespace inside the legend border\nlabelspacing       the vertical space between the legend entries\nhandlelength       the length of the legend handles\nhandletextpad      the pad between the legend handle and text\nborderaxespad      the pad between the axes and legend border\ncolumnspacing      the spacing between columns\n\n'
In [10]:
#function copied from internet to plot with axes in radians
def create_pi_labels(a, b, step):

    max_denominator = int(1/step)
    # i added this line and the .limit_denominator to solve an 
    # issue with floating point precision
    # because of floataing point precision Fraction(1/3) would be 
    # Fraction(6004799503160661, 18014398509481984)

    values = np.arange(a, b+step/10, step)
    fracs = [Fraction(x).limit_denominator(max_denominator) for x in values]
    ticks = values*np.pi

    labels = []

    for frac in fracs:
        if frac.numerator==0:
            labels.append(r"$0$")
        elif frac.numerator<0:
            if frac.denominator==1 and abs(frac.numerator)==1:
                labels.append(r"$-\pi$")
            elif frac.denominator==1:
                labels.append(r"$-{}\pi$".format(abs(frac.numerator)))
            else:
                labels.append(r"$-\frac{{{}}}{{{}}} \pi$".format(abs(frac.numerator), frac.denominator))
        else:
            if frac.denominator==1 and frac.numerator==1:
                labels.append(r"$\pi$")
            elif frac.denominator==1:
                labels.append(r"${}\pi$".format(frac.numerator))
            else:
                labels.append(r"$\frac{{{}}}{{{}}} \pi$".format(frac.numerator, frac.denominator))

    return ticks, labels
In [11]:
#from matplotlib import font_manager
#font_manager.OSXInstalledFonts()

#print  mpl.font_manager.OSXFontDirectories
In [12]:
b=0
minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging
i = 0
for c in [2000]:
    n=Waters.index(c)

    theta2 = array(angles_w[(b, c)])

    end=endpoint(theta2)
    start=best_start(theta2,t,beg,minblocksize)
    print 'start=',start
    slope, intercept, delete1, delete2, delete3 = stats.linregress(t[start:end],theta2[start:end])
    thetamean, errrorbar =blockAverage((theta2[start:end]),blocksNum)
    print "for SAM ",b,"% and ", c, " molecules:"
    if start == None:
        print " Not equilibrated!"
        sampling = 0
        start = beg
        end = endpoint(theta2)
        symb='<'
    else:
        sampling = (end-start)/2.
        print "start point=", start, "length of sampling =", sampling, "ns"
        print "slope=", slope, "intercept=",intercept
        print "average=",thetamean, "error=",errrorbar
        symb='>'
m=str(round(slope*(end-start)/2,3))
print('error '+str(round(2*errrorbar,3))+' deg '+symb + m + ' deg (shift)')
num of good intervals= 52
start= 10
for SAM  0 % and  2000  molecules:
start point= 10 length of sampling = 55.0 ns
slope= 0.0057410771595 intercept= 134.836684409
average= 135.022408718 error= 0.195193442012
error 0.39 deg >0.316 deg (shift)

In [13]:
b=0

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 

errortheta_w[b]=[]
theta_w[b]=[]
start_theta_w[b]=[]
end_theta_w[b]=[]
equil_theta_w[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    theta2 = array(angles_w[(b, c)])

    end=endpoint(theta2)
    start=best_start(theta2,t,beg,minblocksize)
    slope, intercept, delete1, delete2, delete3 = stats.linregress(t[start:end],theta2[start:end])
    thetamean, errrorbar =blockAverage((theta2[start:end]),blocksNum)
    print "for SAM ",b,"% and ", c, " molecules:"
    if start == None:
        print " Not equilibrated!"
        sampling = 0
        start = beg
        end = endpoint(theta2)
        symb='<'
    else:
        sampling = (end-start)/2.
        print "start point=", start, "length of sampling =", sampling, "ns"
        print "slope=", slope, "intercept=",intercept
        print "average=",thetamean, "error=",errrorbar
        symb='>'

    ax = plt.subplot(Nrows, Ncolumns, i+1)

    m=str(round(abs(slope)*(end-start)/2.,3))
    
    line1, = ax.plot(t[start:end], func(t[start:end], slope, intercept),'r', label='error '+str(round(2*errrorbar,3))+'$^\circ $'+symb+m+'$^\circ $ shift')
    ax.plot(t[start:end], theta2[start:end],'-',linewidth=0.5)   
    ax.plot(t[start:end], func(t[start:end], 0, thetamean))
    ax.plot(t[start:end], func(t[start:end], 0, thetamean+errrorbar),'g-')
    ax.plot(t[start:end], func(t[start:end], 0, thetamean-errrorbar),'g-')   
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    plt.ylabel(r'$\theta_{mic}$ [deg]',fontsize=14)
    #ax.set_ylabel(r'R$\ _{BASE} $ [nm]')
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    # Create a legend for the first line.
    #first_legend = plt.legend(handles=[line1],bbox_to_anchor=(0.2, 0.9), loc=3,borderaxespad=0.,borderpad=0.2,fontsize=11)
    first_legend = plt.legend(loc=1,borderaxespad=0.,borderpad=0.2,fontsize=11)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)

    
    # SAVING RESULTS
    if sampling != 0:
        (errortheta_w[b]).append(errrorbar)
        (theta_w[b]).append(thetamean)
        (start_theta_w[b]).append(start)
        (end_theta_w[b]).append(end)
        (equil_theta_w[b]).append(c)
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('equil_theta_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
num of good intervals= 19
for SAM  0 % and  1000  molecules:
start point= 23 length of sampling = 30.5 ns
slope= 0.00342730549085 intercept= 133.21425652
average= 133.305283775 error= 0.146813381649
num of good intervals= 52
for SAM  0 % and  2000  molecules:
start point= 10 length of sampling = 55.0 ns
slope= 0.0057410771595 intercept= 134.836684409
average= 135.022408718 error= 0.195193442012
num of good intervals= 23
for SAM  0 % and  3000  molecules:
start point= 11 length of sampling = 34.0 ns
slope= -0.00180539849189 intercept= 136.235067577
average= 136.190157349 error= 0.080599031913
num of good intervals= 8
for SAM  0 % and  4000  molecules:
start point= 21 length of sampling = 27.5 ns
slope= -0.00548832347994 intercept= 136.888509175
average= 136.757732164 error= 0.0997674132984
num of good intervals= 13
for SAM  0 % and  5000  molecules:
start point= 35 length of sampling = 42.5 ns
slope= -0.00670096077952 intercept= 137.687394482
average= 137.429041522 error= 0.150446740242
for SAM  0 % and  6500  molecules:
 Not equilibrated!
num of good intervals= 21
for SAM  0 % and  7000  molecules:
start point= 76 length of sampling = 42.0 ns
slope= -0.00509756955993 intercept= 137.876329286
average= 137.576847074 error= 0.128212799128
num of good intervals= 16
for SAM  0 % and  8000  molecules:
start point= 92 length of sampling = 34.0 ns
slope= -0.00291560983216 intercept= 137.538917478
average= 137.355080075 error= 0.0578789536767
for SAM  0 % and  9000  molecules:
 Not equilibrated!
for SAM  0 % and  10000  molecules:
 Not equilibrated!

In [14]:
b=0

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 

errortheta_m[b]=[]
theta_m[b]=[]
start_theta_m[b]=[]
end_theta_m[b]=[]
equil_theta_m[b]=[]

errortheta_s[b]=[]
theta_s[b]=[]
start_theta_s[b]=[]
end_theta_s[b]=[]
equil_theta_s[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    print c," molec:"
    n=Waters.index(c)

    theta2_w = array(angles_w[(b, c)])
    theta2_m = array(angles_m[(b, c)])
    theta2_s = array(angles_s[(b, c)])
    
    end=endpoint(theta2_m)
    start=best_start(theta2_m,t,beg,minblocksize)
    thetamean, errrorbar =blockAverage((theta2_m[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (errortheta_m[b]).append(errrorbar)
        (theta_m[b]).append(thetamean)
        (start_theta_m[b]).append(start)
        (end_theta_m[b]).append(end)
        (equil_theta_m[b]).append(c)
        
    end=endpoint(theta2_s)
    start=best_start(theta2_s,t,beg,minblocksize)
    thetamean, errrorbar =blockAverage((theta2_s[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (errortheta_s[b]).append(errrorbar)
        (theta_s[b]).append(thetamean)
        (start_theta_s[b]).append(start)
        (end_theta_s[b]).append(end)
        (equil_theta_s[b]).append(c)
    
    ax = plt.subplot(Nrows, Ncolumns, i+1)
    
    ####line1, = ax.plot(t, theta2_w,'-',linewidth=0.5)
    #ax.plot(t, theta2_w,'-',linewidth=0.5)
    #ax.plot(t, theta2_m,'-',linewidth=0.5)
    #ax.plot(t, theta2_s,'-',linewidth=0.5)
    
    mylines=ax.plot(t, theta2_w,'o', t, theta2_m, 'o', t,theta2_s,'o')
    plt.setp(mylines, markersize=3.0)
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    plt.ylabel(r'$\theta_{mic}$ [deg]',fontsize=14)
    #ax.set_ylabel(r'R$\ _{BASE} $ [nm]')
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec.',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('theta_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
1000  molec:
num of good intervals= 20
num of good intervals= 20
2000  molec:
num of good intervals= 52
num of good intervals= 52
3000  molec:
num of good intervals= 23
num of good intervals= 24
4000  molec:
num of good intervals= 8
num of good intervals= 9
5000  molec:
num of good intervals= 13
num of good intervals= 13
6500  molec:
7000  molec:
num of good intervals= 21
num of good intervals= 21
8000  molec:
num of good intervals= 16
num of good intervals= 16
9000  molec:
10000  molec:

In [15]:
b=5

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 

errortheta_w[b]=[]
theta_w[b]=[]
start_theta_w[b]=[]
end_theta_w[b]=[]
equil_theta_w[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    theta2 = array(angles_w[(b, c)])

    end=endpoint(theta2)
    start=best_start(theta2,t,beg,minblocksize)
    slope, intercept, delete1, delete2, delete3 = stats.linregress(t[start:end],theta2[start:end])
    thetamean, errrorbar =blockAverage((theta2[start:end]),blocksNum)
    print "for SAM ",b,"% and ", c, " molecules:"
    if start == None:
        print " Not equilibrated!"
        sampling = 0
        start = beg
        end = endpoint(theta2)
        symb='<'
    else:
        sampling = (end-start)/2.
        print "start point=", start, "length of sampling =", sampling, "ns"
        print "slope=", slope, "intercept=",intercept
        print "average=",thetamean, "error=",errrorbar
        symb='>'

    ax = plt.subplot(Nrows, Ncolumns, i+1)

    m=str(round(abs(slope)*(end-start)/2.,3))
    line1, = ax.plot(t[start:end], func(t[start:end], slope, intercept),'r', label='error '+str(round(2*errrorbar,3))+'$^\circ $'+symb+m+'$^\circ $ shift')
    ax.plot(t[start:end], theta2[start:end],'-',linewidth=0.5)
    #ax.plot(t, theta2,'-',linewidth=0.5)
    ax.plot(t[start:end], func(t[start:end], 0, thetamean))
    ax.plot(t[start:end], func(t[start:end], 0, thetamean+errrorbar),'g-')
    ax.plot(t[start:end], func(t[start:end], 0, thetamean-errrorbar),'g-')   
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    plt.ylabel(r'$\theta_{mic}$ [deg]',fontsize=14)
    #ax.set_ylabel(r'R$\ _{BASE} $ [nm]')
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    # Create a legend for the first line.
    #first_legend = plt.legend(handles=[line1],bbox_to_anchor=(0.2, 0.9), loc=3,borderaxespad=0.,borderpad=0.2,fontsize=11)
    first_legend = plt.legend(loc=1,borderaxespad=0.,borderpad=0.2,fontsize=11)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)

    
    # SAVING RESULTS
    if sampling != 0:
        (errortheta_w[b]).append(errrorbar)
        (theta_w[b]).append(thetamean)
        (start_theta_w[b]).append(start)
        (end_theta_w[b]).append(end)
        (equil_theta_w[b]).append(c)
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('equil_theta_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
num of good intervals= 49
for SAM  5 % and  1000  molecules:
start point= 21 length of sampling = 49.5 ns
slope= -0.00206583284544 intercept= 127.544613831
average= 127.472309682 error= 0.057870774187
num of good intervals= 39
for SAM  5 % and  2000  molecules:
start point= 30 length of sampling = 45.0 ns
slope= 0.00360507444083 intercept= 126.653894389
average= 126.788183412 error= 0.0843874157528
num of good intervals= 42
for SAM  5 % and  3000  molecules:
start point= 18 length of sampling = 51.0 ns
slope= 0.00484780245488 intercept= 130.623871527
average= 130.789908762 error= 0.167072250734
num of good intervals= 36
for SAM  5 % and  4000  molecules:
start point= 33 length of sampling = 43.5 ns
slope= 0.0029291170748 intercept= 133.56940767
average= 133.680714119 error= 0.0725236149221
num of good intervals= 4
for SAM  5 % and  5000  molecules:
start point= 84 length of sampling = 18.0 ns
slope= 0.00499486249755 intercept= 135.463434078
average= 135.716923349 error= 0.0605684869664
num of good intervals= 27
for SAM  5 % and  6500  molecules:
start point= 77 length of sampling = 41.5 ns
slope= -0.00270234559548 intercept= 129.152169752
average= 128.99073193 error= 0.0578582744324
num of good intervals= 42
for SAM  5 % and  7000  molecules:
start point= 56 length of sampling = 52.0 ns
slope= -0.00127427373672 intercept= 128.920459794
average= 128.853191138 error= 0.0417167632464
for SAM  5 % and  8000  molecules:
 Not equilibrated!
num of good intervals= 20
for SAM  5 % and  9000  molecules:
start point= 71 length of sampling = 24.5 ns
slope= -0.000278743723061 intercept= 130.395573092
average= 130.382113706 error= 0.00790925147783
num of good intervals= 12
for SAM  5 % and  10000  molecules:
start point= 123 length of sampling = 38.0 ns
slope= -0.00451276746715 intercept= 130.631185303
average= 130.270838123 error= 0.0877654170018

In [16]:
b=5

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 

errortheta_m[b]=[]
theta_m[b]=[]
start_theta_m[b]=[]
end_theta_m[b]=[]
equil_theta_m[b]=[]

errortheta_s[b]=[]
theta_s[b]=[]
start_theta_s[b]=[]
end_theta_s[b]=[]
equil_theta_s[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    print c," molec:"
    n=Waters.index(c)

    theta2_w = array(angles_w[(b, c)])
    theta2_m = array(angles_m[(b, c)])
    theta2_s = array(angles_s[(b, c)])
    
    end=endpoint(theta2_m)
    start=best_start(theta2_m,t,beg,minblocksize)
    thetamean, errrorbar =blockAverage((theta2_m[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (errortheta_m[b]).append(errrorbar)
        (theta_m[b]).append(thetamean)
        (start_theta_m[b]).append(start)
        (end_theta_m[b]).append(end)
        (equil_theta_m[b]).append(c)
        
    end=endpoint(theta2_s)
    start=best_start(theta2_s,t,beg,minblocksize)
    thetamean, errrorbar =blockAverage((theta2_s[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (errortheta_s[b]).append(errrorbar)
        (theta_s[b]).append(thetamean)
        (start_theta_s[b]).append(start)
        (end_theta_s[b]).append(end)
        (equil_theta_s[b]).append(c)
    
    ax = plt.subplot(Nrows, Ncolumns, i+1)
    
    ####line1, = ax.plot(t, theta2_w,'-',linewidth=0.5)
    #ax.plot(t, theta2_w,'-',linewidth=0.5)
    #ax.plot(t, theta2_m,'-',linewidth=0.5)
    #ax.plot(t, theta2_s,'-',linewidth=0.5)
    
    mylines=ax.plot(t, theta2_w,'o', t, theta2_m, 'o', t,theta2_s,'o')
    plt.setp(mylines, markersize=3.0)
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    plt.ylabel(r'$\theta_{mic}$ [deg]',fontsize=14)
    #ax.set_ylabel(r'R$\ _{BASE} $ [nm]')
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec.',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('theta_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
1000  molec:
num of good intervals= 50
num of good intervals= 47
2000  molec:
num of good intervals= 39
num of good intervals= 39
3000  molec:
num of good intervals= 42
num of good intervals= 42
4000  molec:
num of good intervals= 38
num of good intervals= 38
5000  molec:
num of good intervals= 4
num of good intervals= 4
6500  molec:
num of good intervals= 26
num of good intervals= 25
7000  molec:
num of good intervals= 42
num of good intervals= 41
8000  molec:
9000  molec:
num of good intervals= 20
num of good intervals= 20
10000  molec:
num of good intervals= 12
num of good intervals= 12

In [17]:
b=11

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
errortheta_w[b]=[]
theta_w[b]=[]
start_theta_w[b]=[]
end_theta_w[b]=[]
equil_theta_w[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    theta2 = array(angles_w[(b, c)])

    end=endpoint(theta2)
    start=best_start(theta2,t,beg,minblocksize)
    slope, intercept, delete1, delete2, delete3 = stats.linregress(t[start:end],theta2[start:end])
    thetamean, errrorbar =blockAverage((theta2[start:end]),blocksNum)
    print "for SAM ",b,"% and ", c, " molecules:"
    if start == None:
        print " Not equilibrated!"
        sampling = 0
        start = beg
        end = endpoint(theta2)
        symb='<'
    else:
        sampling = (end-start)/2.
        print "start point=", start, "length of sampling =", sampling, "ns"
        print "slope=", slope, "intercept=",intercept
        print "average=",thetamean, "error=",errrorbar
        symb='>'

    ax = plt.subplot(Nrows, Ncolumns, i+1)

    m=str(round(abs(slope)*(end-start)/2.,3))
    line1, = ax.plot(t[start:end], func(t[start:end], slope, intercept),'r', label='error '+str(round(2*errrorbar,3))+'$^\circ $'+symb+m+'$^\circ $ shift')
    ax.plot(t[start:end], theta2[start:end],'-',linewidth=0.5)
    #ax.plot(t, theta2,'-',linewidth=0.5)
    ax.plot(t[start:end], func(t[start:end], 0, thetamean))
    ax.plot(t[start:end], func(t[start:end], 0, thetamean+errrorbar),'g-')
    ax.plot(t[start:end], func(t[start:end], 0, thetamean-errrorbar),'g-')   
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    plt.ylabel(r'$\theta_{mic}$ [deg]',fontsize=14)
    #ax.set_ylabel(r'R$\ _{BASE} $ [nm]')
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    # Create a legend for the first line.
    #first_legend = plt.legend(handles=[line1],bbox_to_anchor=(0.2, 0.9), loc=3,borderaxespad=0.,borderpad=0.2,fontsize=11)
    first_legend = plt.legend(loc=1,borderaxespad=0.,borderpad=0.2,fontsize=11)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)

    
    # SAVING RESULTS
    if sampling != 0:
        (errortheta_w[b]).append(errrorbar)
        (theta_w[b]).append(thetamean)
        (start_theta_w[b]).append(start)
        (end_theta_w[b]).append(end)
        (equil_theta_w[b]).append(c)
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('equil_theta_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
num of good intervals= 16
for SAM  11 % and  1000  molecules:
start point= 27 length of sampling = 46.5 ns
slope= -0.0797570993895 intercept= 113.886339734
average= 110.975205607 error= 1.89764360515
num of good intervals= 8
for SAM  11 % and  2000  molecules:
start point= 12 length of sampling = 14.0 ns
slope= 0.0058354143919 intercept= 114.538636235
average= 114.603156761 error= 0.185192368479
for SAM  11 % and  3000  molecules:
 Not equilibrated!
num of good intervals= 20
for SAM  11 % and  4000  molecules:
start point= 32 length of sampling = 44.0 ns
slope= -0.00350566547883 intercept= 112.979232864
average= 112.848411932 error= 0.0801360695192
num of good intervals= 4
for SAM  11 % and  5000  molecules:
start point= 11 length of sampling = 54.5 ns
slope= -0.00147938576794 intercept= 117.416392974
average= 117.368159263 error= 0.0462693380055
num of good intervals= 3
for SAM  11 % and  6500  molecules:
start point= 37 length of sampling = 21.5 ns
slope= -0.00953263947901 intercept= 112.618654386
average= 112.340829558 error= 0.111406990798
num of good intervals= 12
for SAM  11 % and  7000  molecules:
start point= 114 length of sampling = 23.0 ns
slope= -0.021659092225 intercept= 117.216179059
average= 115.735088254 error= 0.270180660303
num of good intervals= 14
for SAM  11 % and  8000  molecules:
start point= 63 length of sampling = 48.5 ns
slope= -0.00302120322651 intercept= 114.350315554
average= 114.182842317 error= 0.092970392899
for SAM  11 % and  9000  molecules:
 Not equilibrated!
num of good intervals= 35
for SAM  11 % and  10000  molecules:
start point= 130 length of sampling = 34.5 ns
slope= -0.00265524977749 intercept= 113.977990223
average= 113.760259741 error= 0.0475681250294

In [18]:
b=11

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 

errortheta_m[b]=[]
theta_m[b]=[]
start_theta_m[b]=[]
end_theta_m[b]=[]
equil_theta_m[b]=[]

errortheta_s[b]=[]
theta_s[b]=[]
start_theta_s[b]=[]
end_theta_s[b]=[]
equil_theta_s[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    print c," molec:"
    n=Waters.index(c)

    theta2_w = array(angles_w[(b, c)])
    theta2_m = array(angles_m[(b, c)])
    theta2_s = array(angles_s[(b, c)])
    
    end=endpoint(theta2_m)
    start=best_start(theta2_m,t,beg,minblocksize)
    thetamean, errrorbar =blockAverage((theta2_m[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (errortheta_m[b]).append(errrorbar)
        (theta_m[b]).append(thetamean)
        (start_theta_m[b]).append(start)
        (end_theta_m[b]).append(end)
        (equil_theta_m[b]).append(c)
        
    end=endpoint(theta2_s)
    start=best_start(theta2_s,t,beg,minblocksize)
    thetamean, errrorbar =blockAverage((theta2_s[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (errortheta_s[b]).append(errrorbar)
        (theta_s[b]).append(thetamean)
        (start_theta_s[b]).append(start)
        (end_theta_s[b]).append(end)
        (equil_theta_s[b]).append(c)
    
    ax = plt.subplot(Nrows, Ncolumns, i+1)
    
    ####line1, = ax.plot(t, theta2_w,'-',linewidth=0.5)
    #ax.plot(t, theta2_w,'-',linewidth=0.5)
    #ax.plot(t, theta2_m,'-',linewidth=0.5)
    #ax.plot(t, theta2_s,'-',linewidth=0.5)
    
    mylines=ax.plot(t, theta2_w,'o', t, theta2_m, 'o', t,theta2_s,'o')
    plt.setp(mylines, markersize=3.0)
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    plt.ylabel(r'$\theta_{mic}$ [deg]',fontsize=14)
    #ax.set_ylabel(r'R$\ _{BASE} $ [nm]')
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec.',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('theta_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
1000  molec:
num of good intervals= 16
num of good intervals= 18
2000  molec:
num of good intervals= 8
num of good intervals= 8
3000  molec:
4000  molec:
num of good intervals= 20
num of good intervals= 20
5000  molec:
num of good intervals= 4
num of good intervals= 4
6500  molec:
num of good intervals= 3
num of good intervals= 3
7000  molec:
num of good intervals= 12
num of good intervals= 12
8000  molec:
num of good intervals= 14
num of good intervals= 14
9000  molec:
10000  molec:
num of good intervals= 35
num of good intervals= 35

In [19]:
b=17

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 

errortheta_w[b]=[]
theta_w[b]=[]
start_theta_w[b]=[]
end_theta_w[b]=[]
equil_theta_w[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    theta2 = array(angles_w[(b, c)])

    end=endpoint(theta2)
    start=best_start(theta2,t,beg,minblocksize)
    slope, intercept, delete1, delete2, delete3 = stats.linregress(t[start:end],theta2[start:end])
    thetamean, errrorbar =blockAverage((theta2[start:end]),blocksNum)
    print "for SAM ",b,"% and ", c, " molecules:"
    if start == None:
        print " Not equilibrated!"
        sampling = 0
        start = beg
        end = endpoint(theta2)
        symb='<'
    else:
        sampling = (end-start)/2.
        print "start point=", start, "length of sampling =", sampling, "ns"
        print "slope=", slope, "intercept=",intercept
        print "average=",thetamean, "error=",errrorbar
        symb='>'

    ax = plt.subplot(Nrows, Ncolumns, i+1)

    m=str(round(abs(slope)*(end-start)/2.,3))
    line1, = ax.plot(t[start:end], func(t[start:end], slope, intercept),'r', label='error '+str(round(2*errrorbar,3))+'$^\circ $'+symb+m+'$^\circ $ shift')
    ax.plot(t[start:end], theta2[start:end],'-',linewidth=0.5)
    #ax.plot(t, theta2,'-',linewidth=0.5)
    ax.plot(t[start:end], func(t[start:end], 0, thetamean))
    ax.plot(t[start:end], func(t[start:end], 0, thetamean+errrorbar),'g-')
    ax.plot(t[start:end], func(t[start:end], 0, thetamean-errrorbar),'g-')   
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    plt.ylabel(r'$\theta_{mic}$ [deg]',fontsize=14)
    #ax.set_ylabel(r'R$\ _{BASE} $ [nm]')
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    # Create a legend for the first line.
    #first_legend = plt.legend(handles=[line1],bbox_to_anchor=(0.2, 0.9), loc=3,borderaxespad=0.,borderpad=0.2,fontsize=11)
    first_legend = plt.legend(loc=1,borderaxespad=0.,borderpad=0.2,fontsize=11)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)

    
    # SAVING RESULTS
    if sampling != 0:
        (errortheta_w[b]).append(errrorbar)
        (theta_w[b]).append(thetamean)
        (start_theta_w[b]).append(start)
        (end_theta_w[b]).append(end)
        (equil_theta_w[b]).append(c)
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('equil_theta_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
for SAM  17 % and  1000  molecules:
 Not equilibrated!
num of good intervals= 36
for SAM  17 % and  2000  molecules:
start point= 25 length of sampling = 47.5 ns
slope= 0.00402294806775 intercept= 100.861819551
average= 101.007614448 error= 0.103658584052
num of good intervals= 13
for SAM  17 % and  3000  molecules:
start point= 25 length of sampling = 47.5 ns
slope= -0.0204712256916 intercept= 102.92954749
average= 102.190565349 error= 0.543688554166
num of good intervals= 20
for SAM  17 % and  4000  molecules:
start point= 33 length of sampling = 43.5 ns
slope= 0.00956714777968 intercept= 100.931655405
average= 101.29520702 error= 0.21430273296
num of good intervals= 24
for SAM  17 % and  5000  molecules:
start point= 50 length of sampling = 35.0 ns
slope= -0.00996213551868 intercept= 102.578089232
average= 102.156464386 error= 0.174598558922
num of good intervals= 40
for SAM  17 % and  6500  molecules:
start point= 34 length of sampling = 43.0 ns
slope= -0.00321358555613 intercept= 105.348424495
average= 105.228734335 error= 0.071421550039
num of good intervals= 5
for SAM  17 % and  7000  molecules:
start point= 135 length of sampling = 12.5 ns
slope= -0.00495097104154 intercept= 103.870671453
average= 103.50932373 error= 0.0481136076704
num of good intervals= 7
for SAM  17 % and  8000  molecules:
start point= 107 length of sampling = 26.5 ns
slope= -0.0156733301421 intercept= 103.842655219
average= 102.814064101 error= 0.218765138489
num of good intervals= 23
for SAM  17 % and  9000  molecules:
start point= 60 length of sampling = 50.0 ns
slope= -0.0151398732873 intercept= 104.893048927
average= 104.062277912 error= 0.407002722401
num of good intervals= 16
for SAM  17 % and  10000  molecules:
start point= 92 length of sampling = 53.5 ns
slope= -0.006605692126 intercept= 102.221355191
average= 101.743039043 error= 0.177207543491

In [20]:
b=17

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 

errortheta_m[b]=[]
theta_m[b]=[]
start_theta_m[b]=[]
end_theta_m[b]=[]
equil_theta_m[b]=[]

errortheta_s[b]=[]
theta_s[b]=[]
start_theta_s[b]=[]
end_theta_s[b]=[]
equil_theta_s[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    print c," molec:"
    n=Waters.index(c)

    theta2_w = array(angles_w[(b, c)])
    theta2_m = array(angles_m[(b, c)])
    theta2_s = array(angles_s[(b, c)])
    
    end=endpoint(theta2_m)
    start=best_start(theta2_m,t,beg,minblocksize)
    thetamean, errrorbar =blockAverage((theta2_m[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (errortheta_m[b]).append(errrorbar)
        (theta_m[b]).append(thetamean)
        (start_theta_m[b]).append(start)
        (end_theta_m[b]).append(end)
        (equil_theta_m[b]).append(c)
        
    end=endpoint(theta2_s)
    start=best_start(theta2_s,t,beg,minblocksize)
    thetamean, errrorbar =blockAverage((theta2_s[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (errortheta_s[b]).append(errrorbar)
        (theta_s[b]).append(thetamean)
        (start_theta_s[b]).append(start)
        (end_theta_s[b]).append(end)
        (equil_theta_s[b]).append(c)
    
    ax = plt.subplot(Nrows, Ncolumns, i+1)
    
    ####line1, = ax.plot(t, theta2_w,'-',linewidth=0.5)
    #ax.plot(t, theta2_w,'-',linewidth=0.5)
    #ax.plot(t, theta2_m,'-',linewidth=0.5)
    #ax.plot(t, theta2_s,'-',linewidth=0.5)
    
    mylines=ax.plot(t, theta2_w,'o', t, theta2_m, 'o', t,theta2_s,'o')
    plt.setp(mylines, markersize=3.0)
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    plt.ylabel(r'$\theta_{mic}$ [deg]',fontsize=14)
    #ax.set_ylabel(r'R$\ _{BASE} $ [nm]')
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec.',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('theta_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
1000  molec:
2000  molec:
num of good intervals= 33
num of good intervals= 31
3000  molec:
num of good intervals= 13
num of good intervals= 13
4000  molec:
num of good intervals= 18
num of good intervals= 17
5000  molec:
num of good intervals= 24
num of good intervals= 24
6500  molec:
num of good intervals= 40
num of good intervals= 40
7000  molec:
num of good intervals= 5
num of good intervals= 5
8000  molec:
num of good intervals= 7
num of good intervals= 7
9000  molec:
num of good intervals= 23
num of good intervals= 23
10000  molec:
num of good intervals= 16
num of good intervals= 15

/Users/burbol2/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/stats/stats.py:3026: RuntimeWarning: invalid value encountered in absolute
  prob = distributions.t.sf(np.abs(t),df)*2

In [21]:
b=33

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 

errortheta_w[b]=[]
theta_w[b]=[]
start_theta_w[b]=[]
end_theta_w[b]=[]
equil_theta_w[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    theta2 = array(angles_w[(b, c)])

    end=endpoint(theta2)
    start=best_start(theta2,t,beg,minblocksize)
    slope, intercept, delete1, delete2, delete3 = stats.linregress(t[start:end],theta2[start:end])
    thetamean, errrorbar =blockAverage((theta2[start:end]),blocksNum)
    print "for SAM ",b,"% and ", c, " molecules:"
    if start == None:
        print " Not equilibrated!"
        sampling = 0
        start = beg
        end = endpoint(theta2)
        symb='<'
    else:
        sampling = (end-start)/2.
        print "start point=", start, "length of sampling =", sampling, "ns"
        print "slope=", slope, "intercept=",intercept
        print "average=",thetamean, "error=",errrorbar
        symb='>'

    ax = plt.subplot(Nrows, Ncolumns, i+1)

    m=str(round(abs(slope)*(end-start)/2.,3))
    line1, = ax.plot(t[start:end], func(t[start:end], slope, intercept),'r', label='error '+str(round(2*errrorbar,3))+'$^\circ $'+symb+m+'$^\circ $ shift')
    ax.plot(t[start:end], theta2[start:end],'-',linewidth=0.5)
    #ax.plot(t, theta2,'-',linewidth=0.5)
    ax.plot(t[start:end], func(t[start:end], 0, thetamean))
    ax.plot(t[start:end], func(t[start:end], 0, thetamean+errrorbar),'g-')
    ax.plot(t[start:end], func(t[start:end], 0, thetamean-errrorbar),'g-')   
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    plt.ylabel(r'$\theta_{mic}$ [deg]',fontsize=14)
    #ax.set_ylabel(r'R$\ _{BASE} $ [nm]')
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    # Create a legend for the first line.
    #first_legend = plt.legend(handles=[line1],bbox_to_anchor=(0.2, 0.9), loc=3,borderaxespad=0.,borderpad=0.2,fontsize=11)
    first_legend = plt.legend(loc=1,borderaxespad=0.,borderpad=0.2,fontsize=11)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)

    
    # SAVING RESULTS
    if sampling != 0:
        (errortheta_w[b]).append(errrorbar)
        (theta_w[b]).append(thetamean)
        (start_theta_w[b]).append(start)
        (end_theta_w[b]).append(end)
        (equil_theta_w[b]).append(c)
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('equil_theta_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
num of good intervals= 36
for SAM  33 % and  1000  molecules:
start point= 14 length of sampling = 70.0 ns
slope= 0.000407447243869 intercept= 67.9973868028
average= 68.0151676871 error= 0.0477011457445
num of good intervals= 24
for SAM  33 % and  2000  molecules:
start point= 10 length of sampling = 69.0 ns
slope= 0.00563153829606 intercept= 68.9528064637
average= 69.1738443418 error= 0.331700455174
num of good intervals= 1
for SAM  33 % and  3000  molecules:
start point= 116 length of sampling = 14.0 ns
slope= -0.0788645500881 intercept= 74.363556675
average= 69.2913619867 error= 0.630536777199
num of good intervals= 20
for SAM  33 % and  4000  molecules:
start point= 13 length of sampling = 65.0 ns
slope= -0.00196465100156 intercept= 69.3510229531
average= 69.2736725838 error= 0.0957857451882
num of good intervals= 53
for SAM  33 % and  5000  molecules:
start point= 10 length of sampling = 67.5 ns
slope= 0.000564885852358 intercept= 69.5299048262
average= 69.5516529315 error= 0.072761071464
num of good intervals= 9
for SAM  33 % and  6500  molecules:
start point= 49 length of sampling = 33.5 ns
slope= 0.0262237835465 intercept= 57.9369796902
average= 59.0016282174 error= 0.44632835029
num of good intervals= 39
for SAM  33 % and  7000  molecules:
start point= 103 length of sampling = 48.0 ns
slope= 0.0126823329169 intercept= 58.9199835961
average= 59.8743291481 error= 0.338474848121
num of good intervals= 44
for SAM  33 % and  8000  molecules:
start point= 54 length of sampling = 53.0 ns
slope= -0.000636603992886 intercept= 69.3878588109
average= 69.354354032 error= 0.0395202573994
num of good intervals= 37
for SAM  33 % and  9000  molecules:
start point= 23 length of sampling = 68.5 ns
slope= -0.00154103421 intercept= 69.8367811846
average= 69.7683036524 error= 0.0837421119346
num of good intervals= 33
for SAM  33 % and  10000  molecules:
start point= 46 length of sampling = 76.5 ns
slope= 0.00490375495162 intercept= 69.4998490569
average= 69.7989781089 error= 0.191664477904

In [87]:
b=33

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 

errortheta_m[b]=[]
theta_m[b]=[]
start_theta_m[b]=[]
end_theta_m[b]=[]
equil_theta_m[b]=[]

errortheta_s[b]=[]
theta_s[b]=[]
start_theta_s[b]=[]
end_theta_s[b]=[]
equil_theta_s[b]=[]
i = 0
for c in Waters:
#for c in [1000]:
    print c," molec:"
    n=Waters.index(c)

    theta2_w = array(angles_w[(b, c)])
    theta2_m = array(angles_m[(b, c)])
    theta2_s = array(angles_s[(b, c)])
    
    end=endpoint(theta2_m)
    start=best_start(theta2_m,t,beg,minblocksize)
    thetamean, errrorbar =blockAverage((theta2_m[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (errortheta_m[b]).append(errrorbar)
        (theta_m[b]).append(thetamean)
        (start_theta_m[b]).append(start)
        (end_theta_m[b]).append(end)
        (equil_theta_m[b]).append(c)
        
    end=endpoint(theta2_s)
    start=best_start(theta2_s,t,beg,minblocksize)
    thetamean, errrorbar =blockAverage((theta2_s[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (errortheta_s[b]).append(errrorbar)
        (theta_s[b]).append(thetamean)
        (start_theta_s[b]).append(start)
        (end_theta_s[b]).append(end)
        (equil_theta_s[b]).append(c)
    
    ax = plt.subplot(Nrows, Ncolumns, i+1)
    
    ####line1, = ax.plot(t, theta2_w,'-',linewidth=0.5)
    #ax.plot(t, theta2_w,'-',linewidth=0.5)
    #ax.plot(t, theta2_m,'-',linewidth=0.5)
    #ax.plot(t, theta2_s,'-',linewidth=0.5)
    
    #print theta2_m
    #print theta2_s
    mylines=ax.plot(t, theta2_w,'ro', t, theta2_m, 'bo', t,theta2_s,'go')
    plt.setp(mylines, markersize=3.0)
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    plt.ylabel(r'$\theta_{mic}$ [deg]',fontsize=14)
    #ax.set_ylabel(r'R$\ _{BASE} $ [nm]')
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec.',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('theta_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
1000  molec:
num of good intervals= 37
num of good intervals= 37
2000  molec:
num of good intervals= 24
num of good intervals= 24
3000  molec:
num of good intervals= 1
num of good intervals= 1
4000  molec:
num of good intervals= 20
num of good intervals= 20
5000  molec:
num of good intervals= 53
num of good intervals= 53
6500  molec:
num of good intervals= 9
num of good intervals= 9
7000  molec:
num of good intervals= 39
num of good intervals= 39
8000  molec:
num of good intervals= 44
num of good intervals= 44
9000  molec:
num of good intervals= 37
num of good intervals= 37
10000  molec:
num of good intervals= 33
num of good intervals= 33

In [23]:
################################ r_base  #################################
b=0

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
error_rbase_w[b]=[]
rbase_w[b]=[]
start_rbase_w[b]=[]
end_rbase_w[b]=[]
equil_rbase_w[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    rbase2 = array(radii_w[(b, c)])

    end=endpoint(rbase2)
    start=best_start(rbase2,t,beg,minblocksize)
    slope, intercept, delete1, delete2, delete3 = stats.linregress(t[start:end],rbase2[start:end])
    rbasemean, errrorbar =blockAverage((rbase2[start:end]),blocksNum)
    print "for SAM ",b,"% and ", c, " molecules:"
    if start == None:
        print " Not equilibrated!"
        sampling = 0
        start = beg
        end = endpoint(rbase2)
        symb='<'
    else:
        sampling = (end-start)/2.
        print "start point=", start, "length of sampling =", sampling, "ns"
        print "slope=", slope, "intercept=",intercept
        print "average=",rbasemean, "error=",errrorbar
        symb='>'

    ax = plt.subplot(Nrows, Ncolumns, i+1)

    m=str(round(abs(slope)*(end-start)/2.,3))
    line1, = ax.plot(t[start:end], func(t[start:end], slope, intercept),'r', label='error '+str(round(2*errrorbar,3))+'$^\circ $'+symb+m+'$^\circ $ shift')
    ax.plot(t[start:end], rbase2[start:end],'-',linewidth=0.5)
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean))
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean+errrorbar),'g-')
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean-errrorbar),'g-')   
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels   
    ax.set_xlabel('time [ns]',fontsize=13)
    #plt.ylabel(r'$\theta_{mic}$ [deg]')
    ax.set_ylabel(r'R$\ _{BASE} $ [nm]',fontsize=14)
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    # Create a legend for the first line.
    #first_legend = plt.legend(handles=[line1],bbox_to_anchor=(0.2, 0.9), loc=3,borderaxespad=0.,borderpad=0.2,fontsize=11)
    first_legend = plt.legend(loc=1,borderaxespad=0.,borderpad=0.2,fontsize=11)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)


    # SAVING RESULTS
    if sampling != 0:        
        (error_rbase_w[b]).append(errrorbar)
        (rbase_w[b]).append(rbasemean)
        (start_rbase_w[b]).append(start)
        (end_rbase_w[b]).append(end)
        (equil_rbase_w[b]).append(c)
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Base Radius $R_{base}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('equil_Rbase_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
num of good intervals= 21
for SAM  0 % and  1000  molecules:
start point= 21 length of sampling = 31.5 ns
slope= -0.000141072210906 intercept= 1.44247039044
average= 1.43880251296 error= 0.00269299288433
num of good intervals= 53
for SAM  0 % and  2000  molecules:
start point= 10 length of sampling = 55.0 ns
slope= -0.00020177371567 intercept= 1.75407577787
average= 1.74754982106 error= 0.006735045407
num of good intervals= 24
for SAM  0 % and  3000  molecules:
start point= 11 length of sampling = 34.0 ns
slope= 8.80992350704e-05 intercept= 1.93562074534
average= 1.93775524314 error= 0.00296534971663
num of good intervals= 9
for SAM  0 % and  4000  molecules:
start point= 21 length of sampling = 27.5 ns
slope= 0.000220842833925 intercept= 2.111586324
average= 2.11684519216 error= 0.00444674095999
num of good intervals= 14
for SAM  0 % and  5000  molecules:
start point= 34 length of sampling = 43.0 ns
slope= 0.000323339434471 intercept= 2.23420139015
average= 2.24658957068 error= 0.00708050482336
for SAM  0 % and  6500  molecules:
 Not equilibrated!
num of good intervals= 21
for SAM  0 % and  7000  molecules:
start point= 76 length of sampling = 42.0 ns
slope= 0.00029341382896 intercept= 2.49006093403
average= 2.50729899648 error= 0.00677196900216
num of good intervals= 17
for SAM  0 % and  8000  molecules:
start point= 92 length of sampling = 34.0 ns
slope= 0.000180986338683 intercept= 2.62094499736
average= 2.6323379207 error= 0.0031660129702
for SAM  0 % and  9000  molecules:
 Not equilibrated!
for SAM  0 % and  10000  molecules:
 Not equilibrated!

In [24]:
b=0

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
error_rbase_m[b]=[]
rbase_m[b]=[]
start_rbase_m[b]=[]
end_rbase_m[b]=[]
equil_rbase_m[b]=[]

error_rbase_s[b]=[]
rbase_s[b]=[]
start_rbase_s[b]=[]
end_rbase_s[b]=[]
equil_rbase_s[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    rbase2_w = array(radii_w[(b, c)])
    rbase2_m = array(radii_m[(b, c)])
    rbase2_s = array(radii_s[(b, c)])
    
    end=endpoint(rbase2_m)
    start=best_start(rbase2_m,t,beg,minblocksize)
    rbasemean, errrorbar =blockAverage((rbase2_m[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (error_rbase_m[b]).append(errrorbar)
        (rbase_m[b]).append(rbasemean)
        (start_rbase_m[b]).append(start)
        (end_rbase_m[b]).append(end)
        (equil_rbase_m[b]).append(c)
        
    end=endpoint(rbase2_s)
    start=best_start(rbase2_s,t,beg,minblocksize)
    rbasemean, errrorbar =blockAverage((rbase2_s[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (error_rbase_s[b]).append(errrorbar)
        (rbase_s[b]).append(rbasemean)
        (start_rbase_s[b]).append(start)
        (end_rbase_s[b]).append(end)
        (equil_rbase_s[b]).append(c)
        
    
    ax = plt.subplot(Nrows, Ncolumns, i+1)
    
    ####line1, = ax.plot(t, rbase2_w,'-',linewidth=0.5)
    #ax.plot(t, rbase2_w,'-',linewidth=0.5)
    #ax.plot(t, rbase2_m,'-',linewidth=0.5)
    #ax.plot(t, rbase2_s,'-',linewidth=0.5)
    
    mylines=ax.plot(t, rbase2_w,'o', t, rbase2_m, 'o', t,rbase2_s,'o')
    plt.setp(mylines, markersize=3.0)
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    ax.set_ylabel(r'R$\ _{BASE} $ [nm]',fontsize=14)
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    ax.set_title(str(c)+' water molec.',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Base Radius $R_{base}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
plt.show()
fig.savefig('rbase_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
num of good intervals= 21
num of good intervals= 20
num of good intervals= 53
num of good intervals= 53
num of good intervals= 24
num of good intervals= 24
num of good intervals= 9
num of good intervals= 9
num of good intervals= 13
num of good intervals= 14
num of good intervals= 21
num of good intervals= 21
num of good intervals= 16
num of good intervals= 16

In [25]:
################################ r_base  #################################
b=5

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold'
error_rbase_w[b]=[]
rbase_w[b]=[]
start_rbase_w[b]=[]
end_rbase_w[b]=[]
equil_rbase_w[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    rbase2 = array(radii_w[(b, c)])

    end=endpoint(rbase2)
    start=best_start(rbase2,t,beg,minblocksize)
    slope, intercept, delete1, delete2, delete3 = stats.linregress(t[start:end],rbase2[start:end])
    rbasemean, errrorbar =blockAverage((rbase2[start:end]),blocksNum)
    print "for SAM ",b,"% and ", c, " molecules:"
    if start == None:
        print " Not equilibrated!"
        sampling = 0
        start = beg
        end = endpoint(rbase2)
        symb='<'
    else:
        sampling = (end-start)/2.
        print "start point=", start, "length of sampling =", sampling, "ns"
        print "slope=", slope, "intercept=",intercept
        print "average=",rbasemean, "error=",errrorbar
        symb='>'

    ax = plt.subplot(Nrows, Ncolumns, i+1)

    m=str(round(abs(slope)*(end-start)/2.,3))
    line1, = ax.plot(t[start:end], func(t[start:end], slope, intercept),'r', label='error '+str(round(2*errrorbar,3))+'$^\circ $'+symb+m+'$^\circ $ shift')
    ax.plot(t[start:end], rbase2[start:end],'-',linewidth=0.5)
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean))
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean+errrorbar),'g-')
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean-errrorbar),'g-')   
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels   
    ax.set_xlabel('time [ns]',fontsize=13)
    #plt.ylabel(r'$\theta_{mic}$ [deg]')
    ax.set_ylabel(r'R$\ _{BASE} $ [nm]',fontsize=14)
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    # Create a legend for the first line.
    #first_legend = plt.legend(handles=[line1],bbox_to_anchor=(0.2, 0.9), loc=3,borderaxespad=0.,borderpad=0.2,fontsize=11)
    first_legend = plt.legend(loc=1,borderaxespad=0.,borderpad=0.2,fontsize=11)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)


    # SAVING RESULTS
    if sampling != 0:        
        (error_rbase_w[b]).append(errrorbar)
        (rbase_w[b]).append(rbasemean)
        (start_rbase_w[b]).append(start)
        (end_rbase_w[b]).append(end)
        (equil_rbase_w[b]).append(c)
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Base Radius $R_{base}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('equil_Rbase_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
num of good intervals= 19
for SAM  5 % and  1000  molecules:
start point= 15 length of sampling = 52.5 ns
slope= -6.2013391385e-05 intercept= 1.57887043843
average= 1.57679298981 error= 0.00164449388221
num of good intervals= 38
for SAM  5 % and  2000  molecules:
start point= 31 length of sampling = 44.5 ns
slope= -5.45839462625e-05 intercept= 2.00473680086
average= 2.00276031901 error= 0.00211987213476
num of good intervals= 45
for SAM  5 % and  3000  molecules:
start point= 18 length of sampling = 51.0 ns
slope= -0.000251990185748 intercept= 2.13962830299
average= 2.13099763913 error= 0.0064816568787
num of good intervals= 38
for SAM  5 % and  4000  molecules:
start point= 35 length of sampling = 42.5 ns
slope= -6.62728286762e-05 intercept= 2.24307261498
average= 2.24052473259 error= 0.0017404945537
num of good intervals= 2
for SAM  5 % and  5000  molecules:
start point= 86 length of sampling = 17.0 ns
slope= -0.000307328239025 intercept= 2.33570038233
average= 2.32009591531 error= 0.00322561807769
num of good intervals= 25
for SAM  5 % and  6500  molecules:
start point= 75 length of sampling = 42.5 ns
slope= 0.000140121418694 intercept= 2.85393003008
average= 2.86217855168 error= 0.00333751982547
num of good intervals= 39
for SAM  5 % and  7000  molecules:
start point= 56 length of sampling = 52.0 ns
slope= 4.59193198059e-05 intercept= 2.93006756945
average= 2.93248130031 error= 0.00192592736653
for SAM  5 % and  8000  molecules:
 Not equilibrated!
num of good intervals= 20
for SAM  5 % and  9000  molecules:
start point= 72 length of sampling = 24.0 ns
slope= 4.02934414546e-05 intercept= 3.11937891868
average= 3.12130293051 error= 0.000702800701977
num of good intervals= 17
for SAM  5 % and  10000  molecules:
start point= 123 length of sampling = 38.0 ns
slope= 0.000218964637716 intercept= 3.22040972044
average= 3.23789082145 error= 0.00469375411219

In [26]:
b=5

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
error_rbase_m[b]=[]
rbase_m[b]=[]
start_rbase_m[b]=[]
end_rbase_m[b]=[]
equil_rbase_m[b]=[]

error_rbase_s[b]=[]
rbase_s[b]=[]
start_rbase_s[b]=[]
end_rbase_s[b]=[]
equil_rbase_s[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    rbase2_w = array(radii_w[(b, c)])
    rbase2_m = array(radii_m[(b, c)])
    rbase2_s = array(radii_s[(b, c)])
    
    end=endpoint(rbase2_m)
    start=best_start(rbase2_m,t,beg,minblocksize)
    rbasemean, errrorbar =blockAverage((rbase2_m[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (error_rbase_m[b]).append(errrorbar)
        (rbase_m[b]).append(rbasemean)
        (start_rbase_m[b]).append(start)
        (end_rbase_m[b]).append(end)
        (equil_rbase_m[b]).append(c)
        
    end=endpoint(rbase2_s)
    start=best_start(rbase2_s,t,beg,minblocksize)
    rbasemean, errrorbar =blockAverage((rbase2_s[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (error_rbase_s[b]).append(errrorbar)
        (rbase_s[b]).append(rbasemean)
        (start_rbase_s[b]).append(start)
        (end_rbase_s[b]).append(end)
        (equil_rbase_s[b]).append(c)
        
    
    ax = plt.subplot(Nrows, Ncolumns, i+1)
    
    ####line1, = ax.plot(t, rbase2_w,'-',linewidth=0.5)
    #ax.plot(t, rbase2_w,'-',linewidth=0.5)
    #ax.plot(t, rbase2_m,'-',linewidth=0.5)
    #ax.plot(t, rbase2_s,'-',linewidth=0.5)
    
    mylines=ax.plot(t, rbase2_w,'o', t, rbase2_m, 'o', t,rbase2_s,'o')
    plt.setp(mylines, markersize=3.0)
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    ax.set_ylabel(r'R$\ _{BASE} $ [nm]',fontsize=14)
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    ax.set_title(str(c)+' water molec.',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Base Radius $R_{base}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
plt.show()
fig.savefig('rbase_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
num of good intervals= 24
num of good intervals= 42
num of good intervals= 38
num of good intervals= 38
num of good intervals= 45
num of good intervals= 46
num of good intervals= 35
num of good intervals= 35
num of good intervals= 2
num of good intervals= 2
num of good intervals= 26
num of good intervals= 25
num of good intervals= 39
num of good intervals= 39
num of good intervals= 20
num of good intervals= 20
num of good intervals= 15
num of good intervals= 15

In [27]:
################################ r_base  #################################
b=11

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
error_rbase_w[b]=[]
rbase_w[b]=[]
start_rbase_w[b]=[]
end_rbase_w[b]=[]
equil_rbase_w[b]=[]

error_rbase_m[b]=[]
rbase_m[b]=[]
start_rbase_m[b]=[]
end_rbase_m[b]=[]
equil_rbase_m[b]=[]

error_rbase_s[b]=[]
rbase_s[b]=[]
start_rbase_s[b]=[]
end_rbase_s[b]=[]
equil_rbase_s[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    rbase2 = array(radii_w[(b, c)])

    end=endpoint(rbase2)
    start=best_start(rbase2,t,beg,minblocksize)
    slope, intercept, delete1, delete2, delete3 = stats.linregress(t[start:end],rbase2[start:end])
    rbasemean, errrorbar =blockAverage((rbase2[start:end]),blocksNum)
    print "for SAM ",b,"% and ", c, " molecules:"
    if start == None:
        print " Not equilibrated!"
        sampling = 0
        start = beg
        end = endpoint(rbase2)
        symb='<'
    else:
        sampling = (end-start)/2.
        print "start point=", start, "length of sampling =", sampling, "ns"
        print "slope=", slope, "intercept=",intercept
        print "average=",rbasemean, "error=",errrorbar
        symb='>'

    ax = plt.subplot(Nrows, Ncolumns, i+1)

    m=str(round(abs(slope)*(end-start)/2.,3))
    line1, = ax.plot(t[start:end], func(t[start:end], slope, intercept),'r', label='error '+str(round(2*errrorbar,3))+'$^\circ $'+symb+m+'$^\circ $ shift')
    ax.plot(t[start:end], rbase2[start:end],'-',linewidth=0.5)
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean))
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean+errrorbar),'g-')
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean-errrorbar),'g-')   
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels   
    ax.set_xlabel('time [ns]',fontsize=13)
    #plt.ylabel(r'$\theta_{mic}$ [deg]')
    ax.set_ylabel(r'R$\ _{BASE} $ [nm]',fontsize=14)
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    # Create a legend for the first line.
    #first_legend = plt.legend(handles=[line1],bbox_to_anchor=(0.2, 0.9), loc=3,borderaxespad=0.,borderpad=0.2,fontsize=11)
    first_legend = plt.legend(loc=1,borderaxespad=0.,borderpad=0.2,fontsize=11)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)


    # SAVING RESULTS
    if sampling != 0:        
        (error_rbase_w[b]).append(errrorbar)
        (rbase_w[b]).append(rbasemean)
        (start_rbase_w[b]).append(start)
        (end_rbase_w[b]).append(end)
        (equil_rbase_w[b]).append(c)
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Base Radius $R_{base}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('equil_Rbase_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
num of good intervals= 26
for SAM  11 % and  1000  molecules:
start point= 12 length of sampling = 54.0 ns
slope= 0.00148972163005 intercept= 1.87326069953
average= 1.92204908292 error= 0.0403338793813
num of good intervals= 8
for SAM  11 % and  2000  molecules:
start point= 12 length of sampling = 14.0 ns
slope= -0.000251331840509 intercept= 2.35594204837
average= 2.35301811435 error= 0.00508274180384
for SAM  11 % and  3000  molecules:
 Not equilibrated!
num of good intervals= 25
for SAM  11 % and  4000  molecules:
start point= 31 length of sampling = 44.5 ns
slope= 0.000108417305686 intercept= 3.01179286334
average= 3.01576064564 error= 0.00269441837268
num of good intervals= 4
for SAM  11 % and  5000  molecules:
start point= 11 length of sampling = 54.5 ns
slope= 4.90773794262e-05 intercept= 3.07512676839
average= 3.07672974992 error= 0.00182618081582
num of good intervals= 3
for SAM  11 % and  6500  molecules:
start point= 37 length of sampling = 21.5 ns
slope= 0.000242344083778 intercept= 3.57001111776
average= 3.57712606997 error= 0.0048421645084
num of good intervals= 10
for SAM  11 % and  7000  molecules:
start point= 114 length of sampling = 23.0 ns
slope= 0.000802378701904 intercept= 3.45820598183
average= 3.51309222918 error= 0.0108841894158
num of good intervals= 12
for SAM  11 % and  8000  molecules:
start point= 62 length of sampling = 49.0 ns
slope= 0.000136430864282 intercept= 3.73857141329
average= 3.74609866126 error= 0.00353744179853
num of good intervals= 1
for SAM  11 % and  9000  molecules:
start point= 59 length of sampling = 10.5 ns
slope= -0.00031720551226 intercept= 3.84706862109
average= 3.83612503092 error= 0.00190349967993
num of good intervals= 35
for SAM  11 % and  10000  molecules:
start point= 129 length of sampling = 35.0 ns
slope= 6.79669563338e-05 intercept= 4.06571527287
average= 4.07123645964 error= 0.00139144702935

In [28]:
b=11

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 

error_rbase_m[b]=[]
rbase_m[b]=[]
start_rbase_m[b]=[]
end_rbase_m[b]=[]
equil_rbase_m[b]=[]

error_rbase_s[b]=[]
rbase_s[b]=[]
start_rbase_s[b]=[]
end_rbase_s[b]=[]
equil_rbase_s[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    rbase2_w = array(radii_w[(b, c)])
    rbase2_m = array(radii_m[(b, c)])
    rbase2_s = array(radii_s[(b, c)])
    
    print 'TEST',rbase2_m[5], rbase2_s[5]
    end=endpoint(rbase2_m)
    start=best_start(rbase2_m,t,beg,minblocksize)
    rbasemean, errrorbar =blockAverage((rbase2_m[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (error_rbase_m[b]).append(errrorbar)
        (rbase_m[b]).append(rbasemean)
        (start_rbase_m[b]).append(start)
        (end_rbase_m[b]).append(end)
        (equil_rbase_m[b]).append(c)
        
    end=endpoint(rbase2_s)
    start=best_start(rbase2_s,t,beg,minblocksize)
    rbasemean, errrorbar =blockAverage((rbase2_s[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (error_rbase_s[b]).append(errrorbar)
        (rbase_s[b]).append(rbasemean)
        (start_rbase_s[b]).append(start)
        (end_rbase_s[b]).append(end)
        (equil_rbase_s[b]).append(c)
        
    
    ax = plt.subplot(Nrows, Ncolumns, i+1)
    
    ####line1, = ax.plot(t, rbase2_w,'-',linewidth=0.5)
    #ax.plot(t, rbase2_w,'-',linewidth=0.5)
    #ax.plot(t, rbase2_m,'-',linewidth=0.5)
    #ax.plot(t, rbase2_s,'-',linewidth=0.5)
    
    mylines=ax.plot(t, rbase2_w,'o', t, rbase2_m, 'o', t,rbase2_s,'o')
    plt.setp(mylines, markersize=3.0)
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    ax.set_ylabel(r'R$\ _{BASE} $ [nm]',fontsize=14)
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    ax.set_title(str(c)+' water molec.',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Base Radius $R_{base}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
plt.show()
fig.savefig('rbase_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
TEST 1.81631012448 1.69942043286
num of good intervals= 23
num of good intervals= 22
TEST 2.30405064842 2.2121632062
num of good intervals= 8
num of good intervals= 8
TEST 2.41569203063 2.30213508381
TEST 2.00581081359 1.80771931317
num of good intervals= 25
num of good intervals= 24
TEST 2.18061866084 1.98606747391
num of good intervals= 4
num of good intervals= 4
TEST 1.97645780423 1.71176897104
num of good intervals= 3
num of good intervals= 3
TEST 1.20945677936 0.617933528302
num of good intervals= 10
num of good intervals= 11
TEST 1.70837213433 1.35075707457
num of good intervals= 12
num of good intervals= 12
TEST 1.21836940412 0.534755746857
num of good intervals= 1
num of good intervals= 1
TEST 1.44701060253 0.927678399654
num of good intervals= 36
num of good intervals= 36

In [29]:
################################ r_base  #################################
b=17

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
error_rbase_w[b]=[]
rbase_w[b]=[]
start_rbase_w[b]=[]
end_rbase_w[b]=[]
equil_rbase_w[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    rbase2 = array(radii_w[(b, c)])

    end=endpoint(rbase2)
    start=best_start(rbase2,t,beg,minblocksize)
    slope, intercept, delete1, delete2, delete3 = stats.linregress(t[start:end],rbase2[start:end])
    rbasemean, errrorbar =blockAverage((rbase2[start:end]),blocksNum)
    print "for SAM ",b,"% and ", c, " molecules:"
    if start == None:
        print " Not equilibrated!"
        sampling = 0
        start = beg
        end = endpoint(rbase2)
        symb='<'
    else:
        sampling = (end-start)/2.
        print "start point=", start, "length of sampling =", sampling, "ns"
        print "slope=", slope, "intercept=",intercept
        print "average=",rbasemean, "error=",errrorbar
        symb='>'

    ax = plt.subplot(Nrows, Ncolumns, i+1)

    m=str(round(abs(slope)*(end-start)/2.,3))
    line1, = ax.plot(t[start:end], func(t[start:end], slope, intercept),'r', label='error '+str(round(2*errrorbar,3))+'$^\circ $'+symb+m+'$^\circ $ shift')
    ax.plot(t[start:end], rbase2[start:end],'-',linewidth=0.5)
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean))
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean+errrorbar),'g-')
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean-errrorbar),'g-')   
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels   
    ax.set_xlabel('time [ns]',fontsize=13)
    #plt.ylabel(r'$\theta_{mic}$ [deg]')
    ax.set_ylabel(r'R$\ _{BASE} $ [nm]',fontsize=14)
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    # Create a legend for the first line.
    #first_legend = plt.legend(handles=[line1],bbox_to_anchor=(0.2, 0.9), loc=3,borderaxespad=0.,borderpad=0.2,fontsize=11)
    first_legend = plt.legend(loc=1,borderaxespad=0.,borderpad=0.2,fontsize=11)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)


    # SAVING RESULTS
    if sampling != 0:        
        (error_rbase_w[b]).append(errrorbar)
        (rbase_w[b]).append(rbasemean)
        (start_rbase_w[b]).append(start)
        (end_rbase_w[b]).append(end)
        (equil_rbase_w[b]).append(c)
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Base Radius $R_{base}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('equil_Rbase_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
for SAM  17 % and  1000  molecules:
 Not equilibrated!
num of good intervals= 13
for SAM  17 % and  2000  molecules:
start point= 52 length of sampling = 34.0 ns
slope= -3.8244163187e-05 intercept= 2.70019123671
average= 2.69844139588 error= 0.00198735592667
num of good intervals= 13
for SAM  17 % and  3000  molecules:
start point= 25 length of sampling = 47.5 ns
slope= 0.000611492317655 intercept= 3.01791724761
average= 3.03998123396 error= 0.0160892153759
num of good intervals= 3
for SAM  17 % and  4000  molecules:
start point= 87 length of sampling = 16.5 ns
slope= -0.000916394365722 intercept= 3.42514247981
average= 3.37794816997 error= 0.00798669476216
num of good intervals= 27
for SAM  17 % and  5000  molecules:
start point= 46 length of sampling = 37.0 ns
slope= 0.000172225586661 intercept= 3.61907292116
average= 3.62622546317 error= 0.00358235881785
num of good intervals= 62
for SAM  17 % and  6500  molecules:
start point= 24 length of sampling = 48.0 ns
slope= 5.69339074677e-05 intercept= 3.8553551537
average= 3.85739054089 error= 0.00173229028182
num of good intervals= 18
for SAM  17 % and  7000  molecules:
start point= 75 length of sampling = 42.5 ns
slope= 0.000175259479675 intercept= 3.98762192889
average= 3.99778560275 error= 0.00398935019087
num of good intervals= 7
for SAM  17 % and  8000  molecules:
start point= 107 length of sampling = 26.5 ns
slope= 0.000629771897617 intercept= 4.18672910198
average= 4.22805288009 error= 0.0088716132392
num of good intervals= 25
for SAM  17 % and  9000  molecules:
start point= 60 length of sampling = 50.0 ns
slope= 0.000676541141029 intercept= 4.32138028116
average= 4.35850029847 error= 0.0177614846558
num of good intervals= 14
for SAM  17 % and  10000  molecules:
start point= 90 length of sampling = 54.5 ns
slope= 0.00029191637686 intercept= 4.59509004314
average= 4.61610266093 error= 0.00830177799524

In [30]:
b=17

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
error_rbase_m[b]=[]
rbase_m[b]=[]
start_rbase_m[b]=[]
end_rbase_m[b]=[]
equil_rbase_m[b]=[]

error_rbase_s[b]=[]
rbase_s[b]=[]
start_rbase_s[b]=[]
end_rbase_s[b]=[]
equil_rbase_s[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    rbase2_w = array(radii_w[(b, c)])
    rbase2_m = array(radii_m[(b, c)])
    rbase2_s = array(radii_s[(b, c)])
    
    end=endpoint(rbase2_m)
    start=best_start(rbase2_m,t,beg,minblocksize)
    rbasemean, errrorbar =blockAverage((rbase2_m[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (error_rbase_m[b]).append(errrorbar)
        (rbase_m[b]).append(rbasemean)
        (start_rbase_m[b]).append(start)
        (end_rbase_m[b]).append(end)
        (equil_rbase_m[b]).append(c)
        
    end=endpoint(rbase2_s)
    start=best_start(rbase2_s,t,beg,minblocksize)
    rbasemean, errrorbar =blockAverage((rbase2_s[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (error_rbase_s[b]).append(errrorbar)
        (rbase_s[b]).append(rbasemean)
        (start_rbase_s[b]).append(start)
        (end_rbase_s[b]).append(end)
        (equil_rbase_s[b]).append(c)
        
    
    ax = plt.subplot(Nrows, Ncolumns, i+1)
    
    ####line1, = ax.plot(t, rbase2_w,'-',linewidth=0.5)
    #ax.plot(t, rbase2_w,'-',linewidth=0.5)
    #ax.plot(t, rbase2_m,'-',linewidth=0.5)
    #ax.plot(t, rbase2_s,'-',linewidth=0.5)
    
    mylines=ax.plot(t, rbase2_w,'o', t, rbase2_m, 'o', t,rbase2_s,'o')
    plt.setp(mylines, markersize=3.0)
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    ax.set_ylabel(r'R$\ _{BASE} $ [nm]',fontsize=14)
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    ax.set_title(str(c)+' water molec.',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Base Radius $R_{base}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
plt.show()
fig.savefig('rbase_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
num of good intervals= 21
num of good intervals= 21
num of good intervals= 13
num of good intervals= 13
num of good intervals= 3
num of good intervals= 3
num of good intervals= 27
num of good intervals= 27
num of good intervals= 62
num of good intervals= 61
num of good intervals= 17
num of good intervals= 15
num of good intervals= 7
num of good intervals= 7
num of good intervals= 26
num of good intervals= 26
num of good intervals= 14
num of good intervals= 14

In [31]:
################################ r_base  #################################
b=33

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 

error_rbase_w[b]=[]
rbase_w[b]=[]
start_rbase_w[b]=[]
end_rbase_w[b]=[]
equil_rbase_w[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    rbase2 = array(radii_w[(b, c)])

    end=endpoint(rbase2)
    start=best_start(rbase2,t,beg,minblocksize)
    slope, intercept, delete1, delete2, delete3 = stats.linregress(t[start:end],rbase2[start:end])
    rbasemean, errrorbar =blockAverage((rbase2[start:end]),blocksNum)
    print "for SAM ",b,"% and ", c, " molecules:"
    if start == None:
        print " Not equilibrated!"
        sampling = 0
        start = beg
        end = endpoint(rbase2)
        symb='<'
    else:
        sampling = (end-start)/2.
        print "start point=", start, "length of sampling =", sampling, "ns"
        print "slope=", slope, "intercept=",intercept
        print "average=",rbasemean, "error=",errrorbar
        symb='>'

    ax = plt.subplot(Nrows, Ncolumns, i+1)

    m=str(round(abs(slope)*(end-start)/2.,3))
    line1, = ax.plot(t[start:end], func(t[start:end], slope, intercept),'r', label='error '+str(round(2*errrorbar,3))+'$^\circ $'+symb+m+'$^\circ $ shift')
    ax.plot(t[start:end], rbase2[start:end],'-',linewidth=0.5)
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean))
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean+errrorbar),'g-')
    ax.plot(t[start:end], func(t[start:end], 0, rbasemean-errrorbar),'g-')   
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels   
    ax.set_xlabel('time [ns]',fontsize=13)
    #plt.ylabel(r'$\theta_{mic}$ [deg]')
    ax.set_ylabel(r'R$\ _{BASE} $ [nm]',fontsize=14)
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    #####ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns')
    ax.set_title(str(c)+' water molec: $\Delta(t)=$'+str(sampling)+'ns',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    # Create a legend for the first line.
    #first_legend = plt.legend(handles=[line1],bbox_to_anchor=(0.2, 0.9), loc=3,borderaxespad=0.,borderpad=0.2,fontsize=11)
    first_legend = plt.legend(loc=1,borderaxespad=0.,borderpad=0.2,fontsize=11)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)


    # SAVING RESULTS
    if sampling != 0:        
        (error_rbase_w[b]).append(errrorbar)
        (rbase_w[b]).append(rbasemean)
        (start_rbase_w[b]).append(start)
        (end_rbase_w[b]).append(end)
        (equil_rbase_w[b]).append(c)
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Base Radius $R_{base}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
###mytitle = plt.suptitle(r'Contact Angle $\theta_{mic}$ vs. Time for SAM'+str(b)+'%', fontsize=17,x=0.7, y=4.65)
plt.show()
fig.savefig('equil_Rbase_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
num of good intervals= 50
for SAM  33 % and  1000  molecules:
start point= 11 length of sampling = 71.5 ns
slope= 1.47557992749e-05 intercept= 2.84920877453
average= 2.84983560404 error= 0.00108706414764
num of good intervals= 26
for SAM  33 % and  2000  molecules:
start point= 10 length of sampling = 69.0 ns
slope= -0.000222841832028 intercept= 3.56495890409
average= 3.55621236218 error= 0.00936471719042
num of good intervals= 9
for SAM  33 % and  3000  molecules:
start point= 21 length of sampling = 61.5 ns
slope= 0.000181407909751 intercept= 4.0170578628
average= 4.0244955871 error= 0.0056607005993
num of good intervals= 20
for SAM  33 % and  4000  molecules:
start point= 13 length of sampling = 65.0 ns
slope= 2.15336036652e-05 intercept= 4.47139997688
average= 4.47227840481 error= 0.00302448754039
num of good intervals= 58
for SAM  33 % and  5000  molecules:
start point= 10 length of sampling = 67.5 ns
slope= -4.35252563292e-05 intercept= 4.80699404172
average= 4.80531831935 error= 0.0029027171928
num of good intervals= 10
for SAM  33 % and  6500  molecules:
start point= 48 length of sampling = 34.0 ns
slope= -0.00121539522384 intercept= 5.77538499679
average= 5.72676454059 error= 0.0209128173274
num of good intervals= 28
for SAM  33 % and  7000  molecules:
start point= 103 length of sampling = 48.0 ns
slope= -0.000691721217151 intercept= 5.87017716
average= 5.81812513841 error= 0.0168034899274
num of good intervals= 49
for SAM  33 % and  8000  molecules:
start point= 54 length of sampling = 53.0 ns
slope= -2.96551039966e-05 intercept= 5.6453925535
average= 5.64380524794 error= 0.00104546975605
num of good intervals= 32
for SAM  33 % and  9000  molecules:
start point= 22 length of sampling = 69.0 ns
slope= 9.41144351496e-05 intercept= 5.86502668755
average= 5.86928536575 error= 0.00395918340833
num of good intervals= 36
for SAM  33 % and  10000  molecules:
start point= 39 length of sampling = 80.0 ns
slope= -0.000247477979724 intercept= 6.09597055904
average= 6.08133773184 error= 0.0100194125237

In [86]:
b=33

minblocksize=20 # The minimum block size is 10ns=20 points
beg = 10 # We leave out the first 5ns=10 points
blocksNum=3 #the function best_start uses also 3 blocks for averaging

Nrows = 5
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
#fig.subplots_adjust(bottom=0.5,top=2.5, left=0.0, right=1.5, hspace = 0.5)
fig.subplots_adjust(bottom=0.0,top=4.5, left=0.0, right=1.5, hspace = 0.35)
matplotlib.rcParams['legend.handlelength'] = 0
matplotlib.rcParams['legend.markerscale'] = 0
matplotlib.rcParams['legend.handletextpad'] = 0
matplotlib.rcParams['legend.markerscale'] = 2 

#matplotlib.rcParams['font.family'] = 'Times New Roman'    
#matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
error_rbase_m[b]=[]
rbase_m[b]=[]
start_rbase_m[b]=[]
end_rbase_m[b]=[]
equil_rbase_m[b]=[]

error_rbase_s[b]=[]
rbase_s[b]=[]
start_rbase_s[b]=[]
end_rbase_s[b]=[]
equil_rbase_s[b]=[]

i = 0
for c in Waters:
#for c in [1000]:
    n=Waters.index(c)

    rbase2_w = array(radii_w[(b, c)])
    rbase2_m = array(radii_m[(b, c)])
    rbase2_s = array(radii_s[(b, c)])
    '''
    if b>=33 and c>=8000:
        print b,c,':'
        print 'radii_m[(b, c)][115:120]=',radii_m[(b, c)][115:120]
        print 'radii_s[(b, c)][115:120]=',radii_s[(b, c)][115:120]
        w=120
        print 'radii_w[z][w]=',radii_w[z][w]
        print 'radii_m[z][w]=',radii_m[z][w]
        print 'radii_s[z][w]=',radii_s[z][w]
        print 'rbase2_w[w]=',rbase2_w[w]
        print 'rbase2_m[w]=',rbase2_m[w]
        print 'rbase2_s[w]=',rbase2_s[w]
     '''
    
    end=endpoint(rbase2_m)
    start=best_start(rbase2_m,t,beg,minblocksize)
    rbasemean, errrorbar =blockAverage((rbase2_m[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (error_rbase_m[b]).append(errrorbar)
        (rbase_m[b]).append(rbasemean)
        (start_rbase_m[b]).append(start)
        (end_rbase_m[b]).append(end)
        (equil_rbase_m[b]).append(c)
        
    end=endpoint(rbase2_s)
    start=best_start(rbase2_s,t,beg,minblocksize)
    rbasemean, errrorbar =blockAverage((rbase2_s[start:end]),blocksNum)
    # SAVING RESULTS
    if start != None:
        (error_rbase_s[b]).append(errrorbar)
        (rbase_s[b]).append(rbasemean)
        (start_rbase_s[b]).append(start)
        (end_rbase_s[b]).append(end)
        (equil_rbase_s[b]).append(c)
        
    
    ax = plt.subplot(Nrows, Ncolumns, i+1)
    
    ####line1, = ax.plot(t, rbase2_w,'-',linewidth=0.5)
    #ax.plot(t, rbase2_w,'-',linewidth=0.5)
    #ax.plot(t, rbase2_m,'-',linewidth=0.5)
    #ax.plot(t, rbase2_s,'-',linewidth=0.5)
    
    mylines1=ax.plot(t, rbase2_w,'r-',label='Water Peak')
    #plt.setp(mylines1, markersize=2.0)
    mylines2=ax.plot(t, rbase2_m, 'b-', label='GDS',markersize=2.5)
    #plt.setp(mylines2, markersize=2.0)
    mylines3=ax.plot(t,rbase2_s,'g-',label='SAM Peak',markersize=2.5)
    #plt.setp(mylines3, markersize=2.0)
    
    matplotlib.rcParams['font.family'] = 'Arial'
    
    # Set common labels
    ax.set_xlabel('time [ns]',fontsize=13)
    ax.set_ylabel(r'R$\ _{BASE} $ [nm]',fontsize=14)
    
    #titles for each subplot:
    matplotlib.rcParams['font.family'] = 'Times New Roman Bold' 
    ax.set_title(str(c)+' water molec.',fontweight='bold',fontsize=15)

    matplotlib.rcParams['font.family'] = 'Arial' 
    
    first_legend = plt.legend(loc=1,borderaxespad=0.,borderpad=0.2,fontsize=11)
    # Add the legend manually to the current Axes.
    rcParams['legend.numpoints'] = 1
    ax = plt.gca().add_artist(first_legend)
    
    i = i+1
 
matplotlib.rcParams['font.family'] = 'Times New Roman' 
mytitle = plt.suptitle(r'Base Radius $R_{base}$ vs. Time for SAM'+str(b)+'%', fontsize=19, fontweight='bold',x=0.7, y=4.65)
plt.show()
fig.savefig('rbase_t_s'+str(b)+'.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
# check for better quality:
# plt.savefig('destination_path.eps', format='eps', dpi=1000)
num of good intervals= 48
num of good intervals= 44
num of good intervals= 26
num of good intervals= 27
num of good intervals= 9
num of good intervals= 8
num of good intervals= 20
num of good intervals= 20
num of good intervals= 57
num of good intervals= 57
num of good intervals= 10
num of good intervals= 10
num of good intervals= 28
num of good intervals= 30
num of good intervals= 50
num of good intervals= 50
num of good intervals= 32
num of good intervals= 32
num of good intervals= 36
num of good intervals= 36

In [33]:
z=33,6500
w=7
print angles_w[z][w]
print angles_m[z][w]
print angles_s[z][w]
print radii_w[z][w]
print radii_m[z][w]
print radii_s[z][w]
60.3038284626
61.1431312889
61.7112275907
5.66037961463
5.70705809286
5.7379588373

In [34]:
for i in SAMs:
    x = rbase_w[i]
    y = theta_w[i]
    print "WaterPeak: SAM",i,"%, len(rbase)=",len(x),"len(theta)=",len(y)
    x = rbase_m[i]
    y = theta_m[i]
    print "MiddlePeak: SAM",i,"%, len(rbase)=",len(x),"len(theta)=",len(y)
    x = rbase_s[i]
    y = theta_s[i]
    print "SAMPeak: SAM",i,"%, len(rbase)=",len(x),"len(theta)=",len(y),"\n"
    #if x.all() != 0:
        #x = 1/x
        #y = cos(radians(theta[i]))  

print equil_theta_w[11], equil_rbase_w[11]
WaterPeak: SAM 0 %, len(rbase)= 7 len(theta)= 7
MiddlePeak: SAM 0 %, len(rbase)= 7 len(theta)= 7
SAMPeak: SAM 0 %, len(rbase)= 7 len(theta)= 7 

WaterPeak: SAM 5 %, len(rbase)= 9 len(theta)= 9
MiddlePeak: SAM 5 %, len(rbase)= 9 len(theta)= 9
SAMPeak: SAM 5 %, len(rbase)= 9 len(theta)= 9 

WaterPeak: SAM 11 %, len(rbase)= 9 len(theta)= 8
MiddlePeak: SAM 11 %, len(rbase)= 9 len(theta)= 8
SAMPeak: SAM 11 %, len(rbase)= 9 len(theta)= 8 

WaterPeak: SAM 17 %, len(rbase)= 9 len(theta)= 9
MiddlePeak: SAM 17 %, len(rbase)= 9 len(theta)= 9
SAMPeak: SAM 17 %, len(rbase)= 9 len(theta)= 9 

WaterPeak: SAM 33 %, len(rbase)= 10 len(theta)= 10
MiddlePeak: SAM 33 %, len(rbase)= 10 len(theta)= 10
SAMPeak: SAM 33 %, len(rbase)= 10 len(theta)= 10 

[1000, 2000, 4000, 5000, 6500, 7000, 8000, 10000] [1000, 2000, 4000, 5000, 6500, 7000, 8000, 9000, 10000]

In [35]:
Nrows = 3
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
fig.subplots_adjust(bottom=0.0,top=1.5, left=0.0, right=2, hspace = 0.25)

k = 0
for i in [0,5,11,17,33]:
#for i in [11]:
    ax = plt.subplot(Nrows, Ncolumns, k+1)
    x = []
    y = []
    x_error = []
    y_error = [] 
    mn=min(len(rbase_w[i]),len(theta_w[i]))
    for j in range(mn):
        x.append(1/rbase_w[i][j])
        x_error.append((rbase_w[i][j]**(-2))*error_rbase_w[i][j])
        y.append(cos(radians(theta_w[i][j])))
        y_error.append(cos(radians(theta_w[i][j]))*radians(errortheta_w[i][j]))
        
    slope, intercept, delete1, delete2, delete3 = stats.linregress(x,y)
    
    yline=np.zeros(len(x))
    m=0
    for l in x:
        yline[m] = func(l, slope, intercept)
        m = m+1

    line2, = ax.plot(x, yline,'b-',label='SAM'+str(i)+'%, y='+str(round(slope,2))+'x+'+str(round(intercept,2)),linewidth=0.5)
    ax.plot(x, y,'x')
    #line2, = ax.errorbar(x, y, fmt='--o')
    ax.errorbar(x, y,xerr=x_error, yerr=y_error, fmt='x')
    ###plt.errorbar(x_WaterPeak, y_WaterPeak,xerr=error_x_WaterPeak, yerr=error_y_WaterPeak,fmt='x', ecolor='red',elinewidth=1.5, color='red',markeredgewidth=1.5)


    for j in range(mn):
        xarrow = x[j]
        yarrow = y[j]
        xtext = ((x[j])*0.01)+x[j]
        ytext = ((y[j])*0.01)+y[j]
        #ax.annotate(str(equil_rbase[i][j]), xy=(xarrow,yarrow), xytext=(xtext,ytext),arrowprops=dict(facecolor='yellow', shrink=0.05,
                                                                                                     # width=1, headwidth=4),
                   # horizontalalignment='right', verticalalignment='top',)
        #ax.annotate(str(equil_rbase[i][j]), xytext=(xtext,ytext), xy=(xarrow,yarrow),horizontalalignment='right', verticalalignment='top',)
    
    # Set common labels
    ax.set_ylabel(r'$cos(\theta_{mic})$ ',fontsize=14)
    ax.set_xlabel('${R \ _{BASE}}^{-1}\ [ nm ^{-1} ]$',fontsize=13)
    
    #titles for each subplot:
    #ax.set_title('ax1 title')
    
    # Create a legend for the first line.
    first_legend = plt.legend(handles=[line2],bbox_to_anchor=(0.2, 0.95), loc=3,borderaxespad=0.)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)
    
    # Set size of subplots
    fig.subplots_adjust(top=4.5) 
    
    plt.grid(True)
    k = k + 1
mytitle = plt.suptitle(r'$cos(\theta_{mic})$ vs 1/R$\ _{BASE} $', fontsize=19, fontweight='bold', x=1.0, y=4.6)
plt.show()
fig.savefig('all_cos_r.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
In [36]:
i=0
x = []
y = []
x_error = []
y_error = []
mn=min(len(rbase_w[i]),len(theta_w[i]))
for j in range(mn):
    x.append(1/rbase_w[i][j])
    x_error.append((rbase_w[i][j]**(-2))*error_rbase_w[i][j])
    y.append(cos(radians(theta_w[i][j])))
    y_error.append(cos(radians(theta_w[i][j]))*radians(errortheta_w[i][j]))
    
slope, intercept, delete1, delete2, delete3 = stats.linregress(x,y)

yline=np.zeros(len(x))
m=0
for l in x:
    yline[m] = func(l, slope, intercept)
    m = m+1

print "Checking vector lengths..."
print len(x),len(y),len(x_error),len(y_error)
print "x=",x
print "y=",y
print "x_error=",x_error
print "y_error=",y_error


# PLOT INSTEAD OF SUBPLOT
plt.figure()
plt.errorbar(x, y, xerr=x_error, yerr=y_error,fmt='x', ecolor='red',elinewidth=1.5, color='black',markeredgewidth=1.5)
plt.plot(x, yline, color='r')
plt.plot(x, y, '-')

plt.xlabel('${R \ _{BASE}}^{-1}\ [ nm ^{-1} ]$',fontsize=13)
plt.ylabel(r' $\cos\Theta_{mic}$ ',fontsize=14)
mytitle = plt.suptitle(r'$\cos\Theta_{mic}$ vs. ${R \ _{BASE}}^{-1}$ for SAM '+str(i)+'%', fontsize=15,fontweight='bold')

# Labels next to each point
for j in range(mn):
    xarrow = x[j]
    yarrow =y[j] + ((y[j])*0.004)
    xtext = x[j]
    ytext = y[j] + ((y[j])*0.01)
    plt.annotate(str(equil_rbase_w[i][j]), xy=(xarrow,yarrow), xytext=(xtext,ytext),arrowprops=dict(facecolor='yellow', shrink=0.05,
                                                                                                  width=1, headwidth=4),
                horizontalalignment='right', verticalalignment='top',)
    #ax.annotate(str(equil_rbase[i][j]), xytext=(xtext,ytext), xy=(xarrow,yarrow),horizontalalignment='right', verticalalignment='top',)
    
plt.grid(True)
#plt.axis([0.3, 0.9, -1.0, 0.0])
plt.show()
fig.savefig('s'+str(i)+'_cos_r.jpg',bbox_extra_artists=(mytitle,), bbox_inches='tight')
Checking vector lengths...
7 7 7 7
x= [0.69502241690159283, 0.57222975159083644, 0.51606104720511969, 0.47240110127287754, 0.44511913215122156, 0.3988355602595039, 0.37989043585087395]
y= [-0.68588546468914235, -0.70738328073894263, -0.72164131648504082, -0.72846342986812451, -0.73644008078516743, -0.73818279826115418, -0.73556618978903443]
x_error= [0.0013008668015999649, 0.0022053696631405221, 0.00078972898435759659, 0.00099234716565788019, 0.0014028677971727585, 0.0010772157827310031, 0.0004569086809495627]
y_error= [-0.001757497067853694, -0.0024098908953235809, -0.0010151461764094244, -0.0012684514059861131, -0.0019337377112172678, -0.0016518578442902671, -0.00074305300995519273]

In [37]:
i=5
x = []
y = []
x_error = []
y_error = []
mn=min(len(rbase_w[i]),len(theta_w[i]))
for j in range(mn):
    x.append(1/rbase_w[i][j])
    x_error.append((rbase_w[i][j]**(-2))*error_rbase_w[i][j])
    y.append(cos(radians(theta_w[i][j])))
    y_error.append(cos(radians(theta_w[i][j]))*radians(errortheta_w[i][j]))
    
slope, intercept, delete1, delete2, delete3 = stats.linregress(x,y)

yline=np.zeros(len(x))
m=0
for l in x:
    yline[m] = func(l, slope, intercept)
    m = m+1

print "Checking vector lengths..."
print len(x),len(y),len(x_error),len(y_error)
print "x=",x
print "y=",y
print "x_error=",x_error
print "y_error=",y_error


# PLOT INSTEAD OF SUBPLOT
plt.figure()
plt.errorbar(x, y, xerr=x_error, yerr=y_error,fmt='x', ecolor='red',elinewidth=1.5, color='black',markeredgewidth=1.5)
plt.plot(x, yline, color='r')
plt.plot(x, y, '-')

plt.xlabel('${R \ _{BASE}}^{-1}\ [ nm ^{-1} ]$',fontsize=13)
plt.ylabel(r' $\cos\Theta_{mic}$ ',fontsize=14)
mytitle = plt.suptitle(r'$\cos\Theta_{mic}$ vs. ${R \ _{BASE}}^{-1}$ for SAM '+str(i)+'%', fontsize=15,fontweight='bold')

# Labels next to each point
for j in range(mn):
    xarrow = x[j]
    yarrow =y[j] + ((y[j])*0.01)
    xtext = x[j]
    ytext = y[j] + ((y[j])*0.03)
    plt.annotate(str(equil_rbase_w[i][j]), xy=(xarrow,yarrow), xytext=(xtext,ytext),arrowprops=dict(facecolor='yellow', shrink=0.05,
                                                                                                  width=1, headwidth=4),
                horizontalalignment='right', verticalalignment='top',)
    #ax.annotate(str(equil_rbase[i][j]), xytext=(xtext,ytext), xy=(xarrow,yarrow),horizontalalignment='right', verticalalignment='top',)
    
plt.grid(True)
#plt.axis([0.3, 0.9, -1.0, 0.0])
plt.show()
fig.savefig('s'+str(i)+'_cos_r.jpg',bbox_extra_artists=(mytitle,), bbox_inches='tight')
Checking vector lengths...
9 9 9 9
x= [0.63419865921533702, 0.49931087135405844, 0.46926377656996171, 0.44632401751889245, 0.43101666332032079, 0.34938421273970477, 0.34100814211343317, 0.32037902832949489, 0.30884302626088017]
y= [-0.60837794039611415, -0.59885844424040835, -0.6532872678469771, -0.69063901966688768, -0.71589899212526587, -0.62919467297646814, -0.6273270478116566, -0.64788213647804715, -0.64640151986632799]
x_error= [0.00066142849563802588, 0.00052850817580044074, 0.0014273158869240604, 0.0003467154414223211, 0.00059924037270044604, 0.00040740880265480634, 0.00022395945475831388, 7.2137376929229234e-05, 0.00044770911203278888]
y_error= [-0.00061448334778279679, -0.00088202162429124791, -0.0019049601060026134, -0.00087419420310105324, -0.00075679079929212426, -0.00063537172144025012, -0.00045675360653099872, -8.9435256644504219e-05, -0.00099015493678179496]

In [38]:
i=11
#we change the label of water 9000
#equil_rbase_w[i].append(10000)
#equil_rbase_w[i][9]=[]
print equil_theta_w[i]
x = []
y = []
x_error = []
y_error = []
mn=min(len(rbase_w[i]),len(theta_w[i]))
for j in range(mn):
    x.append(1/rbase_w[i][j])
    x_error.append((rbase_w[i][j]**(-2))*error_rbase_w[i][j])
    y.append(cos(radians(theta_w[i][j])))
    y_error.append(cos(radians(theta_w[i][j]))*radians(errortheta_w[i][j]))
    
slope, intercept, delete1, delete2, delete3 = stats.linregress(x,y)

yline=np.zeros(len(x))
m=0
for l in x:
    yline[m] = func(l, slope, intercept)
    m = m+1

print "Checking vector lengths..."
print len(x),len(y),len(x_error),len(y_error)
print "x=",x
print "y=",y
print "x_error=",x_error
print "y_error=",y_error


# PLOT INSTEAD OF SUBPLOT
plt.figure()
plt.errorbar(x, y, xerr=x_error, yerr=y_error,fmt='.', ecolor='red',elinewidth=1.5, color='black',markeredgewidth=1.5)
plt.plot(x, yline, color='r')
plt.plot(x, y, '-')

plt.xlabel('${R \ _{BASE}}^{-1}\ [ nm ^{-1} ]$',fontsize=13)
plt.ylabel(r' $\cos\Theta_{mic}$ ',fontsize=14)
mytitle = plt.suptitle(r'$\cos\Theta_{mic}$ vs. ${R \ _{BASE}}^{-1}$ for SAM '+str(i)+'%', fontsize=15,fontweight='bold')

# Labels next to each point
for j in range(mn):
    xarrow = x[j]
    yarrow =y[j] + ((y[j])*0.01)
    xtext = x[j]
    ytext = y[j] + ((y[j])*0.03)
    plt.annotate(str(equil_theta_w[i][j]), xy=(xarrow,yarrow), xytext=(xtext,ytext),arrowprops=dict(facecolor='yellow', shrink=0.05,
                                                                                                  width=1, headwidth=4),
                horizontalalignment='right', verticalalignment='top',)
    #ax.annotate(str(equil_rbase[i][j]), xytext=(xtext,ytext), xy=(xarrow,yarrow),horizontalalignment='right', verticalalignment='top',)
    
plt.grid(True)
#plt.axis([0.3, 0.9, -1.0, 0.0])
plt.show()
fig.savefig('s'+str(i)+'_cos_r.jpg',bbox_extra_artists=(mytitle,), bbox_inches='tight')
[1000, 2000, 4000, 5000, 6500, 7000, 8000, 10000]
Checking vector lengths...
8 8 8 8
x= [0.52027807660505088, 0.42498610355054112, 0.33159130232852829, 0.32502042144751137, 0.2795540275740977, 0.28464951523189935, 0.26694438412479476, 0.26067972027483111]
y= [-0.35796391486259216, -0.41633088678041441, -0.3882943740901818, -0.45970633192370103, -0.38011537703383302, -0.43421082780860987, -0.40964987374110245, -0.40291058481051722]
x_error= [0.010917948648155621, 0.00091801020204501996, 0.00029625882229983297, 0.00019291459004903294, 0.00037841735628603754, 0.00088189521902375903, 0.00025207564126689412, 0.00012935025842694328]
y_error= [-0.011855811015829509, -0.0013456715948167235, -0.00054308336880058846, -0.000371236901492742, -0.00073910348495557999, -0.0020475394377250261, -0.00066471405112887301, -0.0003345045871931614]

In [39]:
i=17
x = []
y = []
x_error = []
y_error = []
mn=min(len(rbase_w[i]),len(theta_w[i]))
for j in range(mn):
    x.append(1/rbase_w[i][j])
    x_error.append((rbase_w[i][j]**(-2))*error_rbase_w[i][j])
    y.append(cos(radians(theta_w[i][j])))
    y_error.append(cos(radians(theta_w[i][j]))*radians(errortheta_w[i][j]))
    
slope, intercept, delete1, delete2, delete3 = stats.linregress(x,y)

yline=np.zeros(len(x))
m=0
for l in x:
    yline[m] = func(l, slope, intercept)
    m = m+1

print "Checking vector lengths..."
print len(x),len(y),len(x_error),len(y_error)
print "x=",x
print "y=",y
print "x_error=",x_error
print "y_error=",y_error


# PLOT INSTEAD OF SUBPLOT
plt.figure()
plt.errorbar(x, y, xerr=x_error, yerr=y_error,fmt='x', ecolor='red',elinewidth=1.5, color='black',markeredgewidth=1.5)
plt.plot(x, yline, color='r')
plt.plot(x, y, '-')

plt.xlabel('${R \ _{BASE}}^{-1}\ [ nm ^{-1} ]$',fontsize=13)
plt.ylabel(r' $\cos\Theta_{mic}$ ',fontsize=14)
mytitle = plt.suptitle(r'$\cos\Theta_{mic}$ vs. ${R \ _{BASE}}^{-1}$ for SAM '+str(i)+'%', fontsize=15,fontweight='bold')

# Labels next to each point
for j in range(mn):
    xarrow = x[j]
    yarrow =y[j] + ((y[j])*0.04)
    xtext = x[j]
    ytext = y[j] + ((y[j])*0.1)
    plt.annotate(str(equil_rbase_w[i][j]), xy=(xarrow,yarrow), xytext=(xtext,ytext),arrowprops=dict(facecolor='yellow', shrink=0.05,
                                                                                                  width=1, headwidth=4),
                horizontalalignment='right', verticalalignment='top',)
    #ax.annotate(str(equil_rbase[i][j]), xytext=(xtext,ytext), xy=(xarrow,yarrow),horizontalalignment='right', verticalalignment='top',)
    
plt.grid(True)
#plt.axis([0.3, 0.9, -1.0, 0.0])
plt.show()
fig.savefig('s'+str(i)+'_cos_r.jpg',bbox_extra_artists=(mytitle,), bbox_inches='tight')
Checking vector lengths...
9 9 9 9
x= [0.37058429415078115, 0.32894939903866188, 0.29603769794024193, 0.27576884288026166, 0.25924261217525213, 0.25013847648830595, 0.23651549031216068, 0.22943671710890917, 0.21663296366077273]
y= [-0.19093944919100952, -0.21116384669508032, -0.19586411192747968, -0.21058205601357233, -0.26267310987152515, -0.23360359430588731, -0.22178785584926505, -0.24297641989646454, -0.20352280381578394]
x_error= [0.00027292899317220618, 0.0017409771053153395, 0.00069994050014106921, 0.00027243285229102868, 0.00011642156865992364, 0.00024961067903061216, 0.00049627429330637821, 0.00093498599319241388, 0.00038960112087279593]
y_error= [-0.00034544451809607942, -0.0020037665509986262, -0.00073258824352416953, -0.00064171085247860158, -0.00032743285491574405, -0.00019616648525158949, -0.00084682417121120768, -0.0017259921274744515, -0.00062946654038258533]

In [40]:
i=33
x = []
y = []
x_error = []
y_error = []
mn=min(len(rbase_w[i]),len(theta_w[i]))
for j in range(mn):
    x.append(1/rbase_w[i][j])
    x_error.append((rbase_w[i][j]**(-2))*error_rbase_w[i][j])
    y.append(cos(radians(theta_w[i][j])))
    y_error.append(cos(radians(theta_w[i][j]))*radians(errortheta_w[i][j]))
    
slope, intercept, delete1, delete2, delete3 = stats.linregress(x,y)

yline=np.zeros(len(x))
m=0
for l in x:
    yline[m] = func(l, slope, intercept)
    m = m+1

print "Checking vector lengths..."
print len(x),len(y),len(x_error),len(y_error)
print "x=",x
print "y=",y
print "x_error=",x_error
print "y_error=",y_error


# PLOT INSTEAD OF SUBPLOT
plt.figure()
plt.errorbar(x, y, xerr=x_error, yerr=y_error,fmt='x', ecolor='red',elinewidth=1.5, color='black',markeredgewidth=1.5)
plt.plot(x, yline, color='r')

plt.xlabel('${R \ _{BASE}}^{-1}\ [ nm ^{-1} ]$',fontsize=13)
plt.ylabel(r' $\cos\Theta_{mic}$ ',fontsize=14)
mytitle = plt.suptitle(r'$\cos\Theta_{mic}$ vs. ${R \ _{BASE}}^{-1}$ for SAM '+str(i)+'%', fontsize=15,fontweight='bold')

# Labels next to each point
for j in range(mn):
    #xarrow = x[j] + ((x[j])*0.02)
    xarrow = x[j]
    yarrow =y[j] + ((y[j])*0.02)
    #xtext = x[j] + ((x[j])*0.1)
    xtext = x[j]
    ytext = y[j] + ((y[j])*0.1)
    plt.annotate(str(equil_rbase_w[i][j]), xy=(xarrow,yarrow), xytext=(xtext,ytext),arrowprops=dict(facecolor='yellow', shrink=0.05,
                                                                                                  width=1, headwidth=4),
                horizontalalignment='right', verticalalignment='top',)
    #ax.annotate(str(equil_rbase[i][j]), xytext=(xtext,ytext), xy=(xarrow,yarrow),horizontalalignment='right', verticalalignment='top',)
    
plt.grid(True)
#plt.axis([0.3, 0.9, -1.0, 0.0])
plt.show()
fig.savefig('s'+str(i)+'_cos_r.jpg',bbox_extra_artists=(mytitle,), bbox_inches='tight')
Checking vector lengths...
10 10 10 10
x= [0.35089743372634596, 0.28119805516503726, 0.24847834426898496, 0.22359967548636037, 0.20810275897288383, 0.17461866869367274, 0.17187667439433041, 0.17718541942348992, 0.1703784937492141, 0.16443750439391175]
y= [0.3743611305456585, 0.35553367551406223, 0.35361586896284297, 0.35390464307024366, 0.34936281623019483, 0.51501371593212475, 0.50189831001964647, 0.35258726949778291, 0.34581732575573731, 0.34531493730852814]
x_error= [0.00013384913121420591, 0.00074049016001417932, 0.00034950007569272041, 0.0001512147436564819, 0.00012570727185909228, 0.00063766692247617352, 0.00049640183018187002, 3.2822180972238577e-05, 0.00011493046658195726, 0.00027092183719099218]
y_error= [0.00031167138314527122, 0.0020582786899797931, 0.0038915224171328072, 0.00059164951153697258, 0.00044366291993272911, 0.0040119049633745886, 0.0029649645349039922, 0.00024320010591938121, 0.00050543815702447836, 0.0011551393092837908]

In [41]:
# IN THE NEXT CELL (BEFORE 3 SEPARATE CELLS) WE REMOVE 2 DATA POINTS FROM SAM33% 
# FROM THE DROPLETS WITH 6500 & 7000 MOLECULES 
i=33
print len(theta_w[i])
print theta_w[i]
print theta_w[i][0:5]
print theta_w[i][7:11]
print len(rbase_w[i])
print rbase_w[i]
print rbase_w[i][0:5]
print rbase_w[i][7:11]
print len(equil_rbase_w[i])
print equil_rbase_w[i]
print len(equil_rbase_m[i])
print equil_rbase_m[i]
print len(equil_rbase_s[i])
print equil_rbase_s[i]
10
[68.015167687090397, 69.173844341775364, 69.291361986659993, 69.27367258378014, 69.551652931477022, 59.001628217399542, 59.874329148119784, 69.354354032005475, 69.768303652399808, 69.798978108901963]
[68.015167687090397, 69.173844341775364, 69.291361986659993, 69.27367258378014, 69.551652931477022]
[69.354354032005475, 69.768303652399808, 69.798978108901963]
10
[2.8498356040411199, 3.5562123621839858, 4.0244955870982109, 4.4722784048092246, 4.8053183193515556, 5.7267645405902394, 5.8181251384102088, 5.6438052479358101, 5.8692853657453625, 6.0813377318381674]
[2.8498356040411199, 3.5562123621839858, 4.0244955870982109, 4.4722784048092246, 4.8053183193515556]
[5.6438052479358101, 5.8692853657453625, 6.0813377318381674]
10
[1000, 2000, 3000, 4000, 5000, 6500, 7000, 8000, 9000, 10000]
10
[1000, 2000, 3000, 4000, 5000, 6500, 7000, 8000, 9000, 10000]
10
[1000, 2000, 3000, 4000, 5000, 6500, 7000, 8000, 9000, 10000]

In [90]:
# RUN ONLY ONCE!!
# IN THE NEXT CELL (BEFORE 3 SEPARATE CELLS) WE REMOVE ONE DATA POINT FROM rbase TO BE ABLE TO PLOT THE RESULTS OF SAM11% 
# THE DATA POINT REMOVED CORRESPONDS TO THE DROPLETS WITH 6500 and 7000 MOLECULES 
i=33
'''
theta_w[i][5:7]=[]
rbase_w[i][5:7]=[]
equil_rbase_w[i][5:7]=[]

theta_m[i][5:7]=[]
rbase_m[i][5:7]=[]
equil_rbase_m[i][5:7]=[]

theta_s[i][5:7]=[]
rbase_s[i][5:7]=[]
equil_rbase_s[i][5:7]=[]

print len(equil_rbase_s[i])

print len(theta_w[i])
print theta_w[i]

print len(rbase_w[i])
print rbase_w[i]

errortheta_w[i][5:7]=[]
error_rbase_w[i][5:7]=[]
errortheta_m[i][5:7]=[]

error_rbase_m[i][5:7]=[]
errortheta_s[i][5:7]=[]
error_rbase_s[i][5:7]=[]

print len(error_rbase_w[i])
print len(error_rbase_m[i])
print len(error_rbase_s[i])

print len(errortheta_w[i])
print len(errortheta_m[i])
print len(errortheta_s[i])
'''
8
6
[68.015167687090397, 69.173844341775364, 69.291361986659993, 69.27367258378014, 69.551652931477022, 69.798978108901963]
6
[2.8498356040411199, 3.5562123621839858, 4.0244955870982109, 4.4722784048092246, 4.8053183193515556, 6.0813377318381674]
6
8
8
6
8
8

In [91]:
# third cell
i=33
x = []
y = []
x_error = []
y_error = []
mn=min(len(rbase_w[i]),len(theta_w[i]))
for j in range(mn):
    x.append(1/rbase_w[i][j])
    x_error.append((rbase_w[i][j]**(-2))*error_rbase_w[i][j])
    y.append(cos(radians(theta_w[i][j])))
    y_error.append(cos(radians(theta_w[i][j]))*radians(errortheta_w[i][j]))

print "Checking vector lengths..."
print len(x),len(y),len(x_error),len(y_error)
print "x=",x
print "y=",y
print "x_error=",x_error
print "y_error=",y_error
    
slope, intercept, delete1, delete2, delete3 = stats.linregress(x,y)

yline=np.zeros(len(x))
m=0
for l in x:
    yline[m] = func(l, slope, intercept)
    m = m+1


# PLOT INSTEAD OF SUBPLOT
plt.figure()
plt.errorbar(x, y, xerr=x_error, yerr=y_error,fmt='x', ecolor='red',elinewidth=1.5, color='black',markeredgewidth=1.5)
plt.plot(x, yline, color='r')

plt.xlabel('${R \ _{BASE}}^{-1}\ [ nm ^{-1} ]$',fontsize=13)
plt.ylabel(r' $\cos\Theta_{mic}$ ',fontsize=14)
mytitle = plt.suptitle(r'$\cos\Theta_{mic}$ vs. ${R \ _{BASE}}^{-1}$ for SAM '+str(i)+'%', fontsize=15,fontweight='bold')

# Labels next to each point
for j in range(mn):
    xarrow = x[j]
    yarrow =y[j] + ((y[j])*0.01)
    xtext = x[j]
    ytext = y[j] + ((y[j])*0.03)
    plt.annotate(str(equil_rbase_w[i][j]), xy=(xarrow,yarrow), xytext=(xtext,ytext),arrowprops=dict(facecolor='yellow', shrink=0.05,
                                                                                                  width=1, headwidth=4),
                horizontalalignment='right', verticalalignment='top',)
    #ax.annotate(str(equil_rbase[i][j]), xytext=(xtext,ytext), xy=(xarrow,yarrow),horizontalalignment='right', verticalalignment='top',)
    
plt.grid(True)
#plt.axis([0.3, 0.9, -1.0, 0.0])
plt.show()
fig.savefig('s'+str(i)+'_cos_r.jpg',bbox_extra_artists=(mytitle,), bbox_inches='tight')
Checking vector lengths...
6 6 6 6
x= [0.35089743372634596, 0.28119805516503726, 0.24847834426898496, 0.22359967548636037, 0.20810275897288383, 0.16443750439391175]
y= [0.3743611305456585, 0.35553367551406223, 0.35361586896284297, 0.35390464307024366, 0.34936281623019483, 0.34531493730852814]
x_error= [0.00013384913121420591, 0.00074049016001417932, 0.00034950007569272041, 0.0001512147436564819, 0.00012570727185909228, 0.00027092183719099218]
y_error= [0.00031167138314527122, 0.0020582786899797931, 0.0038915224171328072, 0.00059164951153697258, 0.00044366291993272911, 0.0011551393092837908]

In [92]:
# WE PLOT THE CORRECTED PLOTS OF SAM 33% AND 11% TOGETHER WITH THE REST.
slopes_all_w=[]
intercepts_all_w=[]

Nrows = 3
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
fig.subplots_adjust(bottom=0.0,top=1.5, left=0.0, right=2, hspace = 0.3)

k = 0
for i in [0,5,11,17,33]:
    ax = plt.subplot(Nrows, Ncolumns, k+1)
    x = []
    y = []
    x_error = []
    y_error = [] 
    #mn=min(len(rbase_w[i]),len(theta_w[i]))
    mn=len(theta_w[i])
    print mn
    for j in range(mn):
        x.append(1/rbase_w[i][j])
        x_error.append((rbase_w[i][j]**(-2))*error_rbase_w[i][j])
        y.append(cos(radians(theta_w[i][j])))
        y_error.append(cos(radians(theta_w[i][j]))*radians(errortheta_w[i][j]))
        
    slope, intercept, delete1, delete2, delete3 = stats.linregress(x,y)
    
    slopes_all_w.append(slope)
    intercepts_all_w.append(intercept)
        
    yline=np.zeros(len(x))
    m=0
    for l in x:
        yline[m] = func(l, slope, intercept)
        m = m+1

    line2, = ax.plot(x, yline,'b-',label='SAM'+str(i)+'%, y='+str(round(slope,2))+'x+'+str(round(intercept,2)),linewidth=0.5)
    #ax.plot(x, y,'x')
    #line2, = ax.errorbar(x, y, fmt='--o')
    ax.errorbar(x, y,xerr=x_error, yerr=y_error, fmt='k.')
    ###plt.errorbar(x_WaterPeak, y_WaterPeak,xerr=error_x_WaterPeak, yerr=error_y_WaterPeak,fmt='x', ecolor='red',elinewidth=1.5, color='red',markeredgewidth=1.5)


    for j in range(mn):
        xarrow = x[j]
        yarrow = y[j]
        xtext = ((x[j])*0.01)+x[j]
        ytext = ((y[j])*0.01)+y[j]
        #ax.annotate(str(equil_rbase[i][j]), xy=(xarrow,yarrow), xytext=(xtext,ytext),arrowprops=dict(facecolor='yellow', shrink=0.05,
                                                                                                     # width=1, headwidth=4),
                   # horizontalalignment='right', verticalalignment='top',)
        #ax.annotate(str(equil_rbase[i][j]), xytext=(xtext,ytext), xy=(xarrow,yarrow),horizontalalignment='right', verticalalignment='top',)
    
    # Set common labels
    ax.set_ylabel(r'$cos(\theta_{mic})$ ',fontsize=15)
    ax.set_xlabel('${R \ _{BASE}}^{-1}\ [ nm ^{-1} ]$',fontsize=13)
    
    #titles for each subplot:
    #ax.set_title('ax1 title')
    
    # Create a legend for the first line.
    first_legend = plt.legend(handles=[line2], loc=1,borderaxespad=0.)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)
    
    # Set size of subplots
    fig.subplots_adjust(top=4.5) 
    
    plt.grid(True)
    k = k + 1
mytitle = plt.suptitle(r'$cos(\theta_{mic})$ vs 1/R$\ _{BASE} $', fontsize=19, fontweight='bold', x=1.0, y=4.6)
plt.show()
fig.savefig('2all_cos_r.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
7
9
8
9
6

In [93]:
# WITHOUT ERRORBARS, AND SYMBOLS AS LEGEND MARKERS (good for white/black print) 
# PLOTS SCALED EQUALLY WITH PRESSURE TENSOR RESULTS AND 3 DEFINITIONS OF Z=0
#(for the 3 def. of z=0)
slopes_all_w=[]
intercepts_all_w=[]
slopes_all_m=[]
intercepts_all_m=[]
slopes_all_s=[]
intercepts_all_s=[]

line1=[]
line2=[]
line3=[]
line4=[]
line5=[]

Nrows = 3
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
fig.subplots_adjust(bottom=0.0,top=1.5, left=0.0, right=2, hspace = 0.3)
matplotlib.rcParams['legend.handlelength'] = 1  #the length of the legend handles
matplotlib.rcParams['legend.markerscale'] = 0.9
matplotlib.rcParams['legend.handletextpad'] = 0.8 #the pad between the legend handle and text
matplotlib.rcParams['font.family'] = 'Times New Roman'
#rcParams['legend.numpoints'] = 1

k = 0
for i in [0,5,11,17,33]:
    ax = plt.subplot(Nrows, Ncolumns, k+1)
    x_w = []
    y_w = []
    x_error_w = []
    y_error_w = []
    #mn=min(len(rbase_w[i]),len(theta_w[i]))
    mn=len(theta_w[i])
    for j in range(mn):
        x_w.append(1/rbase_w[i][j])
        x_error_w.append((rbase_w[i][j]**(-2))*error_rbase_w[i][j])
        y_w.append(cos(radians(theta_w[i][j])))
        y_error_w.append(cos(radians(theta_w[i][j]))*radians(errortheta_w[i][j])) 
    slope, intercept, delete1, delete2, delete3 = stats.linregress(x_w,y_w)
    slopes_all_w.append(slope)
    intercepts_all_w.append(intercept)
    m=0
    l_w=arange(0,1.3,0.1)
    yline_w=np.zeros(len(l_w))
    for l in l_w:
        yline_w[m] = func(l, slope, intercept)
        m = m+1

    x_m = []
    y_m = []
    x_error_m = []
    y_error_m = [] 
    #mn=min(len(rbase_m[i]),len(theta_m[i]))
    mn=len(theta_m[i])
    for j in range(mn):
        x_m.append(1/rbase_m[i][j])
        x_error_m.append((rbase_m[i][j]**(-2))*error_rbase_m[i][j])
        y_m.append(cos(radians(theta_m[i][j])))
        y_error_m.append(cos(radians(theta_m[i][j]))*radians(errortheta_m[i][j]))
    slope, intercept, delete1, delete2, delete3 = stats.linregress(x_m,y_m)
    slopes_all_m.append(slope)
    intercepts_all_m.append(intercept)
    m=0
    l_m=arange(0,1.3,0.1)
    yline_m=np.zeros(len(l_m))
    for l in l_m:
        yline_m[m] = func(l, slope, intercept)
        m = m+1
    
    x_s = []
    y_s = []
    x_error_s = []
    y_error_s = []
    #mn=min(len(rbase_s[i]),len(theta_s[i]))
    mn=len(theta_s[i])
    for j in range(mn):
        x_s.append(1/rbase_s[i][j])
        x_error_s.append((rbase_s[i][j]**(-2))*error_rbase_s[i][j])
        y_s.append(cos(radians(theta_s[i][j])))
        y_error_s.append(cos(radians(theta_s[i][j]))*radians(errortheta_s[i][j]))
    slope, intercept, delete1, delete2, delete3 = stats.linregress(x_s,y_s)
    slopes_all_s.append(slope)
    intercepts_all_s.append(intercept)        
    m=0
    #for l in max(x_s):
    l_s=arange(0,1.3,0.1)
    yline_s=np.zeros(len(l_s))
    for l in l_s:
        yline_s[m] = func(l, slope, intercept)
        m = m+1
    if intercept<0:
        signtext='x'
    else:
        signtext='x+'
    ##############matplotlib.rcParams['legend.markerscale'] = 15
    #plot of the fitted lines
    line2, = ax.plot(x_w, y_w,'c^',markersize=8,label='Water Peak, y = '+str(round(slopes_all_w[k],2))+signtext+str(round(intercepts_all_w[k],2)))
    line3, =ax.plot(x_m, y_m,'o',markersize=7,color='orange',label='Middle Point, y = '+str(round(slopes_all_m[k],2))+signtext+str(round(intercepts_all_m[k],2)))
    line4, =ax.plot(x_s, y_s,'gs',markersize=7,label='SAM Peak, y = '+str(round(slopes_all_s[k],2))+signtext+str(round(intercepts_all_s[k],2)))
    
    #plot of the data points
    ax.plot(l_w, yline_w,'c')
    ax.plot(l_m, yline_m,color='orange')
    ax.plot(l_s, yline_s,'g')
    
    #plot of a black line through y=0
    zeropoints=np.zeros(len(l_w))
    ax.plot(l_w, zeropoints,'k--')

    ##############matplotlib.rcParams['legend.markerscale'] = 0.7
    #plot of pressure tensor point
    #line5, =ax.plot(0,cos(ptensor[k]),'ro',markersize=10, label="Pressure Tensor Method")
    
    ax.set_xlim([0,1.1])
    ax.set_ylim([-1,+1])

    # Set common labels
    ax.set_ylabel(r'$cos(\theta_{mic})$ ',fontsize=15)
    ax.set_xlabel('${R \ _{BASE}}^{-1}\ [ nm ^{-1} ]$',fontsize=13)
    
    #titles for each subplot:
    mysubtitle=ax.set_title('SAM '+str(i)+'%',fontweight='bold',fontsize=14)
    plt.setp(mysubtitle, color='darkgreen')         #set the color of title to white

    # Create a legend for the first line.
    #first_legend = plt.legend(handles=[line2,line3,line4,line5], loc=1,borderaxespad=0.)
    first_legend = plt.legend(handles=[line2,line3,line4], loc=1,borderaxespad=0.)
    first_legend.markerscale=0.7
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)

    
    # Set size of subplots
    fig.subplots_adjust(top=4.5) 
    
    #plt.grid(True)
    k = k + 1
#r'$cos(\theta_{mic})$ vs 1/R$\ _{BASE} $'
titletext='On the same Scale: '+ r'$cos(\theta_{mic})$' +' vs 1/R$\ _{BASE} $ for three Definitions \n of $z=0$ together with Pressure Tensor Results'
mytitle = plt.suptitle(titletext, fontsize=20, fontweight='bold', x=1.0, y=4.7)

plt.show()
fig.savefig('scaled_all_cos_r.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
In [94]:
# WITH ERRORBARS: PLOTS SCALED EQUALLY WITH PRESSURE TENSOR RESULTS AND 3 DEFINITIONS OF Z=0
#(for the 3 def. of z=0)
slopes_all_w=[]
intercepts_all_w=[]
slopes_all_m=[]
intercepts_all_m=[]
slopes_all_s=[]
intercepts_all_s=[]

line1=[]
line2=[]
line3=[]
line4=[]
line5=[]

Nrows = 3
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
fig.subplots_adjust(bottom=0.0,top=1.5, left=0.0, right=2, hspace = 0.3)
matplotlib.rcParams['legend.handlelength'] = 1  #the length of the legend handles
matplotlib.rcParams['legend.markerscale'] = 1.5
matplotlib.rcParams['legend.handletextpad'] = 0.7 #the pad between the legend handle and text
matplotlib.rcParams['font.family'] = 'Times New Roman'
#rcParams['legend.numpoints'] = 1

k = 0
for i in [0,5,11,17,33]:
    ax = plt.subplot(Nrows, Ncolumns, k+1)
    x_w = []
    y_w = []
    x_error_w = []
    y_error_w = []
    #mn=min(len(rbase_w[i]),len(theta_w[i]))
    mn=len(theta_w[i])
    for j in range(mn):
        x_w.append(1/rbase_w[i][j])
        x_error_w.append((rbase_w[i][j]**(-2))*error_rbase_w[i][j])
        y_w.append(cos(radians(theta_w[i][j])))
        y_error_w.append(cos(radians(theta_w[i][j]))*radians(errortheta_w[i][j])) 
    slope, intercept, delete1, delete2, delete3 = stats.linregress(x_w,y_w)
    slopes_all_w.append(slope)
    intercepts_all_w.append(intercept)
    m=0
    l_w=arange(0,1.3,0.1)
    yline_w=np.zeros(len(l_w))
    for l in l_w:
        yline_w[m] = func(l, slope, intercept)
        m = m+1

    x_m = []
    y_m = []
    x_error_m = []
    y_error_m = [] 
    #mn=min(len(rbase_m[i]),len(theta_m[i]))
    mn=len(theta_m[i])
    for j in range(mn):
        x_m.append(1/rbase_m[i][j])
        x_error_m.append((rbase_m[i][j]**(-2))*error_rbase_m[i][j])
        y_m.append(cos(radians(theta_m[i][j])))
        y_error_m.append(cos(radians(theta_m[i][j]))*radians(errortheta_m[i][j]))
    slope, intercept, delete1, delete2, delete3 = stats.linregress(x_m,y_m)
    slopes_all_m.append(slope)
    intercepts_all_m.append(intercept)
    m=0
    l_m=arange(0,1.3,0.1)
    yline_m=np.zeros(len(l_m))
    for l in l_m:
        yline_m[m] = func(l, slope, intercept)
        m = m+1
    
    x_s = []
    y_s = []
    x_error_s = []
    y_error_s = []
    #mn=min(len(rbase_s[i]),len(theta_s[i]))
    mn=len(theta_s[i])
    for j in range(mn):
        x_s.append(1/rbase_s[i][j])
        x_error_s.append((rbase_s[i][j]**(-2))*error_rbase_s[i][j])
        y_s.append(cos(radians(theta_s[i][j])))
        y_error_s.append(cos(radians(theta_s[i][j]))*radians(errortheta_s[i][j]))
    slope, intercept, delete1, delete2, delete3 = stats.linregress(x_s,y_s)
    slopes_all_s.append(slope)
    intercepts_all_s.append(intercept)        
    m=0
    #for l in max(x_s):
    l_s=arange(0,1.3,0.1)
    yline_s=np.zeros(len(l_s))
    for l in l_s:
        yline_s[m] = func(l, slope, intercept)
        m = m+1
    if intercept<0:
        signtext='x'
    else:
        signtext='x+'
    #plot of the fitted lines
    #line2, = ax.plot(x_w, y_w,'c^',markersize=8,label='Water Peak, y = '+str(round(slopes_all_w[k],2))+signtext+str(round(intercepts_all_w[k],2)),linewidth=0.9)
    #line3, =ax.plot(x_m, y_m,'bo',label='Middle Point, y = '+str(round(slopes_all_m[k],2))+signtext+str(round(intercepts_all_m[k],2)),linewidth=0.9)
    #line4, =ax.plot(x_s, y_s,'gs',label='SAM Peak, y = '+str(round(slopes_all_s[k],2))+signtext+str(round(intercepts_all_s[k],2)),linewidth=0.9)

    matplotlib.rcParams['legend.markerscale'] = 1.5
    #plot of the data points
    line2, = ax.plot(l_w, yline_w,'c',label='Water Peak, y = '+str(round(slopes_all_w[k],2))+signtext+str(round(intercepts_all_w[k],2)),linewidth=0.9)
    line3, = ax.plot(l_m, yline_m,color='orange',label='Middle Point, y = '+str(round(slopes_all_m[k],2))+signtext+str(round(intercepts_all_m[k],2)),linewidth=0.9)
    line4, = ax.plot(l_s, yline_s,'g',label='SAM Peak, y = '+str(round(slopes_all_s[k],2))+signtext+str(round(intercepts_all_s[k],2)),linewidth=0.9)
    
    #plot of a black line through y=0
    zeropoints=np.zeros(len(l_w))
    ax.plot(l_w, zeropoints,'k--')
    
    #plots of errorbars
    ax.errorbar(x_w, y_w,xerr=x_error_w, yerr=y_error_w, fmt='c.',markersize=8,label='Water Peak, y = '+str(round(slopes_all_w[k],2))+signtext+str(round(intercepts_all_w[k],2)),linewidth=0.9)
    ax.errorbar(x_m, y_m,xerr=x_error_m, yerr=y_error_m, fmt='.',color='orange' ,label='Middle Point, y = '+str(round(slopes_all_m[k],2))+signtext+str(round(intercepts_all_m[k],2)),linewidth=0.9)
    ax.errorbar(x_s, y_s,xerr=x_error_s, yerr=y_error_s, fmt='g.',label='SAM Peak, y = '+str(round(slopes_all_s[k],2))+signtext+str(round(intercepts_all_s[k],2)),linewidth=0.9)
    #other options for errorbars: ecolor='red',elinewidth=1.5, color='red',markeredgewidth=1.5
    
    ax.set_xlim([0,1.1])
    ax.set_ylim([-1,+1])

    # Set common labels
    ax.set_ylabel(r'$cos(\theta_{mic})$ ',fontsize=15)
    ax.set_xlabel('${R \ _{BASE}}^{-1}\ [ nm ^{-1} ]$',fontsize=13)
    
    #titles for each subplot:
    mysubtitle=ax.set_title('SAM '+str(i)+'%',fontweight='bold',fontsize=14)
    plt.setp(mysubtitle, color='darkgreen')         #set the color of title to white

    
    matplotlib.rcParams['legend.markerscale'] = 0.8
    #plot of pressure tensor point
    #line5, =ax.plot(0,cos(ptensor[k]),'ro',markersize=10, label="Pressure Tensor")
    
    # Create a legend for the first line.
    #first_legend = plt.legend(handles=[line2,line3,line4,line5], loc=1,borderaxespad=0.)
    first_legend = plt.legend(handles=[line2,line3,line4], loc=1,borderaxespad=0.)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)
    
    
    # Set size of subplots
    fig.subplots_adjust(top=4.5) 
    
    #plt.grid(True)
    k = k + 1
#r'$cos(\theta_{mic})$ vs 1/R$\ _{BASE} $'
titletext='On the Same Scale: '+ r'$cos(\theta_{mic})$' +' vs 1/R$\ _{BASE} $ for Three Definitions \n OF $z=0$ Together with Pressure Tensor Results'
mytitle = plt.suptitle(titletext, fontsize=20, fontweight='bold', x=1.0, y=4.7)

plt.show()
fig.savefig('scaled_errors_all_cos_r.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
In [118]:
# WITH ERRORBARS: PLOTS SCALED EQUALLY WITH PRESSURE TENSOR RESULTS AND 3 DEFINITIONS OF Z=0
#(for the 3 def. of z=0)
slopes_all_w=[]
intercepts_all_w=[]
slopes_all_m=[]
intercepts_all_m=[]
slopes_all_s=[]
intercepts_all_s=[]

line1=[]
line2=[]
line3=[]
line4=[]
line5=[]

Nrows = 3
Ncolumns = 2
fig, ax = plt.subplots(Ncolumns, Nrows)
fig.subplots_adjust(bottom=0.0,top=1.5, left=0.0, right=2, hspace = 0.3)
matplotlib.rcParams['legend.handlelength'] = 1  #the length of the legend handles
matplotlib.rcParams['legend.markerscale'] = 1.5
matplotlib.rcParams['legend.handletextpad'] = 0.7 #the pad between the legend handle and text
matplotlib.rcParams['font.family'] = 'Times New Roman'
#rcParams['legend.numpoints'] = 1

k = 0
for i in [0,5,11,17,33]:
    ax = plt.subplot(Nrows, Ncolumns, k+1)
    x_w = []
    y_w = []
    x_error_w = []
    y_error_w = []
    #mn=min(len(rbase_w[i]),len(theta_w[i]))
    mn=len(theta_w[i])
    for j in range(mn):
        x_w.append(1/rbase_w[i][j])
        x_error_w.append((rbase_w[i][j]**(-2))*error_rbase_w[i][j])
        y_w.append(cos(radians(theta_w[i][j])))
        y_error_w.append(cos(radians(theta_w[i][j]))*radians(errortheta_w[i][j])) 
    slope, intercept, delete1, delete2, delete3 = stats.linregress(x_w,y_w)
    slopes_all_w.append(slope)
    intercepts_all_w.append(intercept)
    m=0
    l_w=arange(0,1.3,0.1)
    yline_w=np.zeros(len(l_w))
    for l in l_w:
        yline_w[m] = func(l, slope, intercept)
        m = m+1

    x_m = []
    y_m = []
    x_error_m = []
    y_error_m = [] 
    #mn=min(len(rbase_m[i]),len(theta_m[i]))
    mn=len(theta_m[i])
    for j in range(mn):
        x_m.append(1/rbase_m[i][j])
        x_error_m.append((rbase_m[i][j]**(-2))*error_rbase_m[i][j])
        y_m.append(cos(radians(theta_m[i][j])))
        y_error_m.append(cos(radians(theta_m[i][j]))*radians(errortheta_m[i][j]))
    slope, intercept, delete1, delete2, delete3 = stats.linregress(x_m,y_m)
    slopes_all_m.append(slope)
    intercepts_all_m.append(intercept)
    m=0
    l_m=arange(0,1.3,0.1)
    yline_m=np.zeros(len(l_m))
    for l in l_m:
        yline_m[m] = func(l, slope, intercept)
        m = m+1
    
    x_s = []
    y_s = []
    x_error_s = []
    y_error_s = []
    #mn=min(len(rbase_s[i]),len(theta_s[i]))
    mn=len(theta_s[i])
    for j in range(mn):
        x_s.append(1/rbase_s[i][j])
        x_error_s.append((rbase_s[i][j]**(-2))*error_rbase_s[i][j])
        y_s.append(cos(radians(theta_s[i][j])))
        y_error_s.append(cos(radians(theta_s[i][j]))*radians(errortheta_s[i][j]))
    slope, intercept, delete1, delete2, delete3 = stats.linregress(x_s,y_s)
    slopes_all_s.append(slope)
    intercepts_all_s.append(intercept)        
    m=0
    #for l in max(x_s):
    l_s=arange(0,1.3,0.1)
    yline_s=np.zeros(len(l_s))
    for l in l_s:
        yline_s[m] = func(l, slope, intercept)
        m = m+1
    if intercept<0:
        signtext='x'
    else:
        signtext='x+'
    #plot of the fitted lines
    #line2, = ax.plot(x_w, y_w,'c^',markersize=8,label='Water Peak, y = '+str(round(slopes_all_w[k],2))+signtext+str(round(intercepts_all_w[k],2)),linewidth=0.9)
    #line3, =ax.plot(x_m, y_m,'bo',label='Middle Point, y = '+str(round(slopes_all_m[k],2))+signtext+str(round(intercepts_all_m[k],2)),linewidth=0.9)
    #line4, =ax.plot(x_s, y_s,'gs',label='SAM Peak, y = '+str(round(slopes_all_s[k],2))+signtext+str(round(intercepts_all_s[k],2)),linewidth=0.9)

    matplotlib.rcParams['legend.markerscale'] = 1.5
    #plot of the data points
    line2, = ax.plot(l_w, yline_w,'c',label='Water Peak, y = '+str(round(slopes_all_w[k],2))+signtext+str(round(intercepts_all_w[k],2)),linewidth=0.9)
    line3, = ax.plot(l_m, yline_m,color='orange',label='Middle Point, y = '+str(round(slopes_all_m[k],2))+signtext+str(round(intercepts_all_m[k],2)),linewidth=0.9)
    line4, = ax.plot(l_s, yline_s,'g',label='SAM Peak, y = '+str(round(slopes_all_s[k],2))+signtext+str(round(intercepts_all_s[k],2)),linewidth=0.9)
    
    #plot of a black line through y=0
    zeropoints=np.zeros(len(l_w))
    ax.plot(l_w, zeropoints,'k--')
    
    #plots of errorbars
    ax.errorbar(x_w, y_w,xerr=x_error_w, yerr=y_error_w, fmt='c.',markersize=8,label='Water Peak, y = '+str(round(slopes_all_w[k],2))+signtext+str(round(intercepts_all_w[k],2)),linewidth=0.9)
    ax.errorbar(x_m, y_m,xerr=x_error_m, yerr=y_error_m, fmt='.',color='orange' ,label='Middle Point, y = '+str(round(slopes_all_m[k],2))+signtext+str(round(intercepts_all_m[k],2)),linewidth=0.9)
    ax.errorbar(x_s, y_s,xerr=x_error_s, yerr=y_error_s, fmt='g.',label='SAM Peak, y = '+str(round(slopes_all_s[k],2))+signtext+str(round(intercepts_all_s[k],2)),linewidth=0.9)
    #other options for errorbars: ecolor='red',elinewidth=1.5, color='red',markeredgewidth=1.5
    
    ax.set_xlim([-0.1,1.1])
    if i==0:
        ax.set_ylim([-0.95,-0.55])
    elif i==5:
        ax.set_ylim([-0.85,-0.55])
    elif i==11:
        ax.set_ylim([-0.6,-0.3])
    elif i==17:
        ax.set_ylim([-0.45,-0.15])
    elif i==33:
        ax.set_ylim([0.22,0.52])

    # Set common labels
    ax.set_ylabel(r'$cos(\theta_{mic})$ ',fontsize=16)
    ax.set_xlabel('${R \ _{BASE}}^{-1}\ [ nm ^{-1} ]$',fontsize=15)
    
    #titles for each subplot:
    mysubtitle=ax.set_title('SAM '+str(i)+'%',fontweight='bold',fontsize=15)
    plt.setp(mysubtitle, color='darkgreen')         #set the color of title to white

    
    matplotlib.rcParams['legend.markerscale'] = 0.8
    #plot of pressure tensor point
    #line5, =ax.plot(0,cos(ptensor[k]),'ro',markersize=10, label="Pressure Tensor")
    
    # Create a legend for the first line.
    #first_legend = plt.legend(handles=[line2,line3,line4,line5], loc=1,borderaxespad=0.)
    first_legend = plt.legend(handles=[line2,line3,line4], loc=1,borderaxespad=0.)
    # Add the legend manually to the current Axes.
    ax = plt.gca().add_artist(first_legend)
    
    
    # Set size of subplots
    fig.subplots_adjust(top=4.5) 
    
    #plt.grid(True)
    k = k + 1
#r'$cos(\theta_{mic})$ vs 1/R$\ _{BASE} $'
titletext='On the Same Scale: '+ r'$cos(\theta_{mic})$' +' vs 1/R$\ _{BASE} $ for Three Definitions \n OF $z=0$ Together with Pressure Tensor Results'
mytitle = plt.suptitle(titletext, fontsize=20, fontweight='bold', x=1.0, y=4.7)

plt.show()
#fig.savefig('scaled_errors_all_cos_r.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
In [53]:
print slopes_all_w
print intercepts_all_w
print len(slopes_all_w), len(intercepts_all_w)

print slopes_all_m
print intercepts_all_m
print len(slopes_all_m), len(intercepts_all_m)

print slopes_all_s
print intercepts_all_s
print len(slopes_all_s), len(intercepts_all_s)
[0.17171391833889232, 0.079573734302370455, 0.14503068179629192, 0.24343566586400942, 0.1311039031148947]
[-0.8072930495786429, -0.68002506881397673, -0.45498147920368509, -0.28586332267783765, 0.32391596774896642]
5 5
[0.11830045572063613, 0.011241331828832227, 0.0093108861910405783, -0.030031788365810084, 0.092959766030251925]
[-0.82399678663367615, -0.69945190930668655, -0.46346144777911824, -0.27607104009266831, 0.31240881626577238]
5 5
[-0.05825809628508917, -0.14115511839436445, -0.16766344393382715, -0.12055496910057809, -0.016802343653108023]
[-0.80738475992576697, -0.68753059597852706, -0.45117050260484926, -0.27630520177933582, 0.32350515628566778]
5 5

In [54]:
# WE CALCULATE THE MACROSCOPIC CONTACT ANGLES (theta_mac_w in radians) AND THE LINE TENSIONS (sigma_w)

theta_mac_w=np.zeros(len(intercepts_all_w))
sigma_w=np.zeros(len(slopes_all_w))
for i in range(len(intercepts_all_w)):
    theta_mac_w[i]=arccos(intercepts_all_w[i])
    sigma_w[i]=(-52.7)*slopes_all_w[i]
print theta_mac_w
print sigma_w

theta_mac_m=np.zeros(len(intercepts_all_m))
sigma_m=np.zeros(len(slopes_all_m))
for i in range(len(intercepts_all_m)):
    theta_mac_m[i]=arccos(intercepts_all_m[i])
    sigma_m[i]=(-52.7)*slopes_all_m[i]
print theta_mac_m
print sigma_m

theta_mac_s=np.zeros(len(intercepts_all_s))
sigma_s=np.zeros(len(slopes_all_s))
for i in range(len(intercepts_all_s)):
    theta_mac_s[i]=arccos(intercepts_all_s[i])
    sigma_s[i]=(-52.7)*slopes_all_s[i]
print theta_mac_s
print sigma_s
[ 2.51034706  2.31859315  2.04314774  1.86070355  1.24093063]
[  9.0493235    4.1935358    7.64311693  12.82905959   6.90917569]
[ 2.53922564  2.34542663  2.05269386  1.8505002   1.25306861]
[ 6.23443402  0.59241819  0.4906837  -1.58267525  4.89897967]
[ 2.51050248  2.32887923  2.03887281  1.85074384  1.24136482]
[-3.07020167 -7.43887474 -8.8358635  -6.35324687 -0.88548351]

In [55]:
x=sigma_w
y=theta_mac_w
print "x=",x
print "y=",y

fig = plt.figure()
ax  = fig.add_subplot(111)

ax.plot(x,y,'o-')

y_pi   = y/np.pi
y_tick,y_label=create_pi_labels(0.4, 0.8, 0.1) # DON'T CHANGE!!!!!!!!!!
#print y_tick
#print y_label
ax.set_yticks(y_tick)

ax.set_yticklabels(y_label, fontsize=20)

ax2 = ax.twinx()
ax2.set_yticks(ax.get_yticks())
ax2.set_ylim((1.199,2.6)) # DON'T CHANGE!!!!!!!!!!

ax.set_ylabel(r'$\theta_{mac}$ [rad]',fontsize=15)
ax.set_xlabel('$\sigma$ [pN]',fontsize=15)
mytitle = plt.suptitle(r'Macroscopic Contact Angle $\theta_{mac}$ vs Line Tension $\sigma $', fontsize=19, fontweight='bold')
#plt.grid(which='major', axis='both') #(use grid if numbers are changed)
plt.show()
fig.savefig('theta_mac_sigma_w.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
x= [  9.0493235    4.1935358    7.64311693  12.82905959   6.90917569]
y= [ 2.51034706  2.31859315  2.04314774  1.86070355  1.24093063]

In [56]:
#plots without complete wetting results
# Two subplots, unpack the axes array immediately
fig, axarr = plt.subplots(2, sharex=True)

# we mark the zero point with a line
zerox=range(36)
linezero=np.zeros(36)
axarr[1].plot(zerox,linezero,'k-')

# we mark the 90 deg. point with a line for the contact angle
zerox=range(36)
lineninety=36*[90]
axarr[0].plot(zerox,lineninety,'k-')

axarr[0].plot(SAMs,degrees(theta_mac_w),'c^--')
axarr[0].plot(SAMs,degrees(theta_mac_m),'o--', color='orange')
axarr[0].plot(SAMs,degrees(theta_mac_s),'gs--')
mytitle1=axarr[0].set_title('Macroscopic Cotact Angle',fontsize=16, fontweight='bold')

axarr[1].plot(SAMs,sigma_w,'c^--',label="Water Peak")
axarr[1].plot(SAMs,sigma_m,'o--', color='orange',label="Middle Point")
axarr[1].plot(SAMs,sigma_s,'gs--',label="SAM Peak")
axarr[1].set_title('Line Tension',fontsize=16, fontweight='bold')

axarr[0].set_ylabel(r'$\theta_{mac}$ [deg]',fontsize=15)
axarr[1].set_ylabel('$\sigma$ [pN]',fontsize=15)
axarr[1].set_xlabel('% OH- coverage percentage',fontsize=13)

first_legend=axarr[1].legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

plt.show()
fig.savefig('subplt_alltheta_mac_sigma.jpg',bbox_extra_artists=(first_legend,mytitle1,), bbox_inches='tight')
In [57]:
# we add the results of SAM 50% and SAM 60% to sigma (in sigma2) and to theta_mac (in theta_mac2)
# and change all the types to list (incl. percentages-> percentages2)

sigma2_w=list(sigma_w)
sigma2_w=sigma2_w+[0,0]
print sigma2_w

theta_mac2_w=list(theta_mac_w)
theta_mac2_w=theta_mac2_w+[0,0]
print theta_mac2_w

sigma2_m=list(sigma_m)
sigma2_m=sigma2_m+[0,0]
print sigma2_m

theta_mac2_m=list(theta_mac_m)
theta_mac2_m=theta_mac2_m+[0,0]
print theta_mac2_m

sigma2_s=list(sigma_s)
sigma2_s=sigma2_s+[0,0]
print sigma2_s

theta_mac2_s=list(theta_mac_s)
theta_mac2_s=theta_mac2_s+[0,0]
print theta_mac2_s

percentages2=list(percentages)
print percentages2
[9.0493234964596265, 4.1935357977349232, 7.6431169306645845, 12.829059591033296, 6.9091756941549507, 0, 0]
[2.5103470648111683, 2.3185931523950183, 2.043147742635155, 1.860703551668679, 1.2409306312736113, 0, 0]
[6.2344340164775245, 0.59241818737945839, 0.49068370226783853, -1.5826752468781915, 4.8989796697942767, 0, 0]
[2.5392256381622214, 2.3454266313120038, 2.0526938625698543, 1.850500197823435, 1.2530686125395492, 0, 0]
[-3.0702016742241995, -7.4388747393830075, -8.835863495312692, -6.3532468716004651, -0.88548351051879293, 0, 0]
[2.5105024828874414, 2.3288792255824657, 2.0388728106330287, 1.8507438361662882, 1.2413648211955797, 0, 0]
[0, 5, 11, 17, 33, 50, 66]

In [60]:
# plots with complete wetting results AND with pressure tensor results
# Two subplots, unpack the axes array immediately
fig, axarr = plt.subplots(2, sharex=True)

# we mark the zero point with a line for the line tension
zerox=range(71)
linezero=np.zeros(71)
axarr[1].plot(zerox,linezero,'k-')

# we mark the 90 deg. point with a line for the contact angle
zerox=range(71)
lineninety=71*[90]
axarr[0].plot(zerox,lineninety,'k-')

axarr[0].plot(percentages2,degrees(theta_mac2_w),'c^--',label='Droplet Method: equal for \n the three z=0 defintions')
#axarr[0].plot(percentages2,degrees(theta_mac2_m),'o--', color='orange',label='Middle Point: Droplet Method')
#axarr[0].plot(percentages2,degrees(theta_mac2_s),'gs--',label='SAM Peak: Droplet Method')
mytitle1=axarr[0].set_title('Macroscopic Cotact Angle',fontsize=16, fontweight='bold')

#line2, =axarr[0].plot(percentages2,ptensor,'rs--',label='Pressure Tensor Method')

axarr[1].plot(percentages2,sigma2_w,'c^--',label="Water Peak")
axarr[1].plot(percentages2,sigma2_m,'o--',color='orange',label="Middle Point")
axarr[1].plot(percentages2,sigma2_s,'gs--',label="SAM Peak")
axarr[1].set_title('Line Tension',fontsize=16, fontweight='bold')

axarr[0].set_ylabel(r'$\theta_{mac}$ [deg]',fontsize=15)
axarr[1].set_ylabel('$\sigma$ [pN]',fontsize=15)
axarr[1].set_xlabel('% OH- coverage percentage',fontsize=13)

first_legend=axarr[1].legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
second_legend=axarr[0].legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

plt.show()
fig.savefig('subplt2_theta_mac_sigma.jpg',bbox_extra_artists=(first_legend,second_legend,mytitle1,), bbox_inches='tight')
In [61]:
x=sigma_w
y=theta_mac_w
print "x=",x
print "y=",y

fig = plt.figure()
ax1 = fig.add_subplot(121)

ax1.plot(x_w,y_w,'c^-')
ax1.plot(x_m,y_m,'o-',color='orange')
ax1.plot(x_s,y_s,'gs-')

ax2  = fig.add_subplot(122)

ax2.plot(x2_w,y2_w,'c^-',label="Water Peak")
ax2.plot(x2_m,y2_m,'o-',color='orange',label="Middle Point")
ax2.plot(x2_s,y2_s,'gs-',label="SAM Peak")

#y_pi   = y/np.pi
#y_tick,y_label=create_pi_labels(0.4, 0.8, 0.1) # DON'T CHANGE!!!!!!!!!!
#print y_tick
#print y_label
#ax.set_yticks(y_tick)

#ax.set_yticklabels(y_label, fontsize=20)

ax1.set_ylabel(r'$\theta_{mac}$ [rad]',fontsize=15)
ax1.set_xlabel('$\sigma$ [pN]',fontsize=15)
ax2.set_xlabel('$\sigma$ [pN]',fontsize=15)

first_legend=ax2.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

mytitle = plt.suptitle(r'Macroscopic Contact Angle $\theta_{mac}$ vs Line Tension $\sigma $', fontsize=19, fontweight='bold')
#plt.grid(which='major', axis='both') #(use grid if numbers are changed)

fig.savefig('theta_mac_sigma2.jpg',bbox_extra_artists=(first_legend,mytitle,), bbox_inches='tight')
plt.show()
x= [  9.0493235    4.1935358    7.64311693  12.82905959   6.90917569]
y= [ 2.51034706  2.31859315  2.04314774  1.86070355  1.24093063]

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-61-de4d3cfd51f0> in <module>()
     13 ax2  = fig.add_subplot(122)
     14 
---> 15 ax2.plot(x2_w,y2_w,'c^-',label="Water Peak")
     16 ax2.plot(x2_m,y2_m,'o-',color='orange',label="Middle Point")
     17 ax2.plot(x2_s,y2_s,'gs-',label="SAM Peak")

NameError: name 'x2_w' is not defined
In []:
# mycell.show()
i=98
#      %save file_name.py _
#     %save?
#    %save file_name2 %history 103
In []: